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SUMMARY 


The objective of this research program is to extend the recent advances 
in robust control system design of multivariable systems to sensor failure 
detection, isolation, and accommodation (DIA), and estimator design. This 
effort provides analysis tools to quantify the trade-off between performance 
robustness and DIA sensitivity, which are to be used to achieve higher levels 
of performance robustness for given levels of DIA sensitivity. An 
innovations-based DIA scheme is used. Estimators, which depend upon a model 
of the process and process Inputs and outputs, are used to generate these 
innovations. Thresholds used to determine failure detection are computed 
based on bounds on modeling errors, noise properties, and the class of 
failures. The applicability of the newly developed tools are demonstrated on 
a multivariable aircraft turbojet engine example. 

A new concept called the threshold selector was developed under this 
program. It represents a significant and 'nnovative tool ror the analysts ana 
synthesis of OIA algorithms. Analytical results were obtained for the SISO 
case to compute optimal thresholds and to determine the size of minimum 
detectable failures, and a computer-aided technique was developed for the 
multivariable case. 

The estimators were made robust by introduction of an internal model and 
by frequency shaping. The internal model provides asymptotically unbiased 
filter estimates. The Incorporation of frequency shaping of the LQG cost 
functional modifies the estimator design to make it suitable for sensor 
failure OIA. 


The results are compared with previous studies which used thresholds that 
were selected empirically. Comparison of these two techniques on a nonlinear 
dynamic engine simulation shows improved performance of the new method 
compared to previous techniques. 
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I. INTRODUCTION 


The use of analytical redundancy (as opposed to hardware redundancy) for 
actuator/sensor failure detection, isolation, and accommodation (DIA) has been 
an active area of research during the last decade. Failure detection is the 
process of determining if a malfunction has occurred In a system. Failures in 
a system are detectable If the outputs following the failure are statistically 
different from the outputs prior to the failure. Failure detection implies 
that a no-failure condition can be differentiated from a failure condition. 

It does not Imply that In the case of a failure, various failures can be dif- 
ferentiated from each other. Failure isolation is the process of differenti- 
ating between various failures which may occur In a system. Any two failures 
may be differentiated from each other if the outputs following them are sta- 
tistically different. In general, more complex data processing Is required 
for isolation than for detection. Failure accommodation refers to the substi- 
tution of a synthesized value for the faulty sensor's output. 

The various techniques that have been developed For DIA can generally be 
thought of as belonging to three categories: failure-sensitive filters, 

multiple-hypothesis filter detectors, and innovations-based detection sys- 
tems. A survey of various techniques up to 1975 can be found in Ref. 1 and 
references to the more recent work is containedlin Refs. 2, 3, and 27. In 
Ref. 1, the issue of robustness of various techniques was pointed out as an 
open area for research. However, little has been done In the way of robust- 
ness analysis for various failure detection schemes. This study seems to be 
the first to address this Important Issue directly. It presents a systematic 
solution based on recent robustness analysis and design techniques [ 14 ] devel- 
oped for multivariable systems. For example, Clark [2] points out the impor- 
tance of robustness but his schemes have no guaranteed robustness properties, 
lelnlnger [4] addresses the problem of parameter uncertainty (D.C. mismatch), 
but does not come up with a remedy to guarantee robustness. The use of adap- 
tivity for failure detection was discussed in a recent report [5]. However, 
these techniques are known to have undesirable robustness properties unless 
high-frequency unmodeled dynamics are taken into account. 



This report addresses the important issue of robustness of sensor failure 
detection, isolation, and accommodation (OIA) techniques. The approach is 
based on extensions of robust control and estimation techniques and unifies 
various OIA methods. It is Important to note that the method has been made 
practical so as to make an immediate impact on applied technology. 

1.1 PROBLEM STATEMENT 

The overall problem addressed In this report Is the design of a robust, 
multivariable control system for a jet engine. A typical block diagram of 
such a system is shown in Figure 1.1. The blocks in the forward loop — the 
actuators, plant dynamics, and sensors — constitute the system. The blocks 
in the feedback system — the control law, DIA logic, and the state estimator 
contain the subsystems which are to be designed. The design objective, 
ideally. Is to select the proper configuration of these feedback subsystems 
such that the closed-loop system exhibits performance robustness with respect 
to system uncertainties and sensor failures. 

In this study, the control law obtained Prom Refs. 5 and 7 was used and 
only the DIA subsystem was designed. In general, a comprehensive design tech- 
nique would also consider the design of the control law as well. 

Performance robustness requirements can be stated as tolerances on: 

(PI) asymptotic behavior, e.g., steady-state command following and 
disturbance rejection; and 

(P2) transient behavior, e.g., speed of response, damping, over- 
shoot, etc. 

(P3) detection and isolation sensitivity, e.g. ability and speed of 
correctly detecting a sensor failure. 

Uncertainties in the system are caused by: 

(Ul) uncertain parameter values in the models of the actuators, 
plant dynamics, and sensors; 

(U2) unmodeled dynamics, e.g., the effect of neglecting high- 
frequency phenomena, neglecting nonlinearities, and 
intentional reduced-order modeling; 
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Figure 1.1 Generalized Block Diagram of a System with Sensor OIA 


3 






(U3) sensor failures of a known type, e.g., slow drift in a sensor 
bias; 

(1)4) DIA reconf iguration of the estimator and/or control law when 
the DIA presumes a sensor failure; and 

(U5) uncertain external signals, e.g., sensor noise and 
disturbances . 

Model uncertainties of the type in (Ul) affect the steady-state 
(asymptotic) behavior and are the predominant cause of estimator bias and 
steady-state regulation errors. Those of the type in (U2) typically have a 
greater effect at higher frequencies and show up In the transient response as 
(possibly) undesirable behavior. A primary goal of this effort is to quantify 
the effect of these uncertainties (U1-U5) on DIA system performance and 
provide a design approach for improving DIA performance (P1-P3). 

1.2 CONTRIBUTIONS 


The general contribution of this research has been the extension of 
recent advances in robust control system design to sensor DIA and estimator 
design. The soecific contributions are: 

- analysis tools with which to quanti fy the trade-off between 
performance robustness and DIA sensitivity; 

- design methods which allow higher levels of performance 
robustness to be achieved for given levels of DIA sensitivity; 

- demonstration of the applicability of these tools using an 
aircraft turbine jet engine multivariable control example. 

The requirement for these goals is explained with the aid of Figure 1.2. 
Plotted (conceptually) are levels of performance robustness against DIA 
sensitivity. A trade-off Is indicated. As a design becomes more robust, it 
becomes less sensitive. Alternately, for a given level of DIA sensitivity, 
there Is a maximum level of performance robustness achievable within the 
estimator and DIA logic design being applied. There are three curves in 
Figure 1.2. The one labeled "current" refers to the current state of the art 
DIA algorithms. The curve labeled "robust" refers to the idea of making the 
0 I A robust. This is the subject of the present study and will result in 
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higher levels of performance robustness. The curve labeled "adaptive" points 
to the fact that even further Improvements in performance may be achievable by 
an adaptive scheme. This requires basic research and would be the subject of 
future studies. 

Curves such as those indicated In Figure 1.2, and more specifically in 
Figure 1.3, constitute a powerful design aid. The quantities in these figures 
will be defined precisely in the body of the report and these figures are 
presented here only to provide a flavor of the results to be discussed. 

Figure 1.3(a) shows the threshold, J th , in an innovations-based DIA scheme 
as a function of a moving detection window, x. This threshold is a fixed 
level against which some measure of the innovations signal is being compared. 

A failure is declared if the measure exceeds the threshold. Figure 1.3(b) 
shows the minimum detectable level of failure which is possible for a given 
DIA technique. These curves are functions of model error bounds, R.M.S. 
noise, and the class of failure signals. Figure 1.3(b) shows that the effect 
of different estimator speeds can also be evaluated. The ability to generate 
these curves is a powerful synthesis tool. 

Plots such as shown in Figure 1.3 provide the means by which to evaluate 
design modifications made in search of the proper balance of robustness and 
DIA sensitivity. Furthermore, they provide information applicable to 
optimizing the search. The sensitivities of the system performance to design 
changes are calculated and the effects of critical parameters are identified. 

This research has developed the analysis tools required to construct 
these curves. In particular: 

- a "threshold selector" has been created which quantifies the 
effects of uncertainty on DIA performance; 

- measures have been derived to quantify the uncertainty and the 
performance robustness. 

This study has also provided advances to robust estimator design to 
achieve higher levels of performance robustness. Specifically, the 
accomplishments are: 
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J: Threshold 

: Threshold for large detection window 
n: RMS noise 

t: length of moving detection window 
Figure 1.3a Threshold vs. Detection Window 




f(-r): smallest size of detectable failure 
t f : time fai lure occurs 

t: length of detection window 

Figure 1.3b Minimum Detectable Bias Shift vs. Detection Window 
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- the development of estimators using the “internal model 

principle 11 to achieve asymptotic convergence despite model error; 

the Incorporation of frequency weighting in an LQG cost 
functional to modify an estimator design to be suitable for 
sensor failure 0 1 A . 

1.3 BACKGROUND 

1.3.1 Previous Programs 

Under subcontract to Pratt & Whitney Aircraft (PWA), Systems Control 
Technology, Inc. (SCT) conducted two previous studies In the area of sensor 
DIA under NASA contracts NAS3-22481 , "Sensor Failure Detection System" and 
NAS3-23282, "Sensor Failure Detection for Jet Engines." The present program 
is an extension of these studies. The ultimate objective of all these 
programs is to provide sensor fail operational control capability while 
minimizing the required sensor hardware redundancy. 

"Sensor Failure Detection System," NA S3-22481 T 8 1) 

The objective of this study was to develop an advanced concept for the 
DIA of sensor failures in gas turbine engine control systems. Five concepts 
were formulated from advanced techniques for sensor DIA. These concepts were 
evaluated by application to a turbofan engine and multivariable control system 
simulation. A simplified version of the simulation was used in the 
preliminary screening process to select one of the DIA concepts. This 
simplified model was also used for the filter designs in the various DIA 
concepts . 

A functional diagram of the selected advanced concept is shown in Figure 
1.4. A normal mode filter, l.e., a filter designed to use all sensor inputs 
with no failures assumed on those Inputs, was used to generate the filter 
residuals and the estimated measurements. The DIA concept used innovations 
from the filter to detect hard failures and used the weighted sum-squared 
residuals (WSSR) technique to detect soft failures. Isolation of soft 
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Figure 1.4 Detection, Isolation, and Accommodation (DIA) Concept [3] 
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failures was accomplished by likelihood ratio (LR) based testing of 
innovations from a bank of Kalman filters, each designed with the assumption 
of one failed input. Accommodation was accomplished by reconfiguring the 
normal mode filter to eliminate the failed sensor from the inputs to the 
filter. The DIA concept is summarized in Table 1.1. The DIA concept selected 
was compared against a baseline DIA concept based on the conventional 
techniques of parameter synthesis. 


Table 1.1 

Advanced Concept for Detecting, Isolating, and Accommodating 

Sensor Failures 

Detection - Innovations testing based on WSSR technique for soft 

failure. Innovations testing against thresholds for 
hard failures. 

Isolation - On-line isolation of hard failures using Innovations 

testing; off-line isolation of soft failures using LR 
technique. Both structures employ bank of Kalman 
filters . 

Accommodation - Reconfiguration and reinitialization of normal mode 

filter. 


The configuration of the multivariable control and the components of the 
DIA logic used in conjunction with the engine simulation is shown in Figure 

1.5. The form of the control law is given by 

u = u b+ c p ( z p - z pb ) + Cj J (Zj - z Ib ) dt (1.1) 

where u is the input vector [WF AH CIVV RCVV 8LC]\ Zp is the estimate of 

the output vector, [N1 N2 PT4 PT6 FTIT] T , and Zj is a subset of the vector 

Zp. ‘ denotes the estimates. Up, z pb> and z^p are the base point vectors 
and Cp and Cj are proportional and integral control gain matrices. The 
proportional part of the control law provides regulation and the Integral part 
provides trim for the fan speed (Nl) and augmentor pressure (PT6). Note that 
the control law uses the estimates for both the proportional and integral 
portions . 
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Figure 1.5 Detailed Block Diagram of DIA Concept [S] 
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As a result of this program, the advanced DIA concept was shown to be a 
viable DIA technique for application to gas turbine engines. While the 
performance of the advanced concept was generally good and its feasibility was 
demonstrated, several problem areas were identified. These Included: 

- steady-state and dynamic mismatch of the simplified nonlinear 
model ; 

- steady-state estimate errors with no failures induced; 

- Instabilities when accommodating failures; 

- accommodation inaccuracies; 

- missed detections and false alarms; and 
limited coverage on the flight envelope. 

The models used In the no-failure filter and the Isolation filters 
contained modeling errors which contributed to the above problems. The hard 
and soft failure-detection thresholds and the soft failure-isolation threshold 
were chosen empirically. They encompassed an estimate of model errors, and 
estimates of sensor noise and bias errors, and a built-in safety factor. 

Since the model errors were large, the thresholds were large and contributed 
to missed detections and false alarms and the overall poor detection 
performance for very slow drift failures. 

"Sensor Failure Detection for Jet Engines, 11 N AS3-23?82 [91 

The objective of this program was to develop refinements to the sensor 
failure OIA algorithm (Figure 1.4) developed under NASA Contract NAS3-22481 . 
These refinements included: 

- improvement of the steady-state accuracy of the simplified model 
of the engine; 

improvement of the dynamic characteristics of the simplified 
model to be comparable with the nonlinear thermodynamic 
simulation; 

- refinements to the DIA algorithm to be compatible with the 
improved model ; 
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- elimination of the steady-state errors with no sensor failures 
caused by biased filter estimates; 

- accommodation inaccuracies in case of failure detection and 
isolation; 

- missed failure detections and false alarms. 

Three revisions were developed and evaluated to address the above 
Improvements. A detailed, nonlinear thermodynamic simulation was used to 
evaluate each revision. As a result, one revision was chosen for detailed 
evaluation and full envelope operation. This DIA algorithm has been 
programmed on a real-time, microprocessor-based controls computer [10, 11] at 
NASA Lewis Research Center in preparation for testing as part of a closed-loop 
control of an engine in a test cell. 


This DIA concept, referred to as Revision 2 in the previous contract [9], 
operated as follows. In the case of no sensor failures, the outputs of the 
normal mode filter were fed to the proportional control part and the sensed 
outputs (N1 and PT6 only) were fed to the integral control. This ensured that 
in steady state, the engine outputs were at the reference point, i.e., there 
were no steady-state hang-off errors. In the case of a failure, the 
synthesized value of the outputs were fed to both the proportional and 
integral controls. Note that if the failed sensor is N1 or PT6, the integral 
logic ensured that the estimates of these measurements were driven to the 
reference point. A detailed block diagram of the overall system is shown in 
Figure 1.6. 


The form of the control law, as shown in Figure 1.6, is given by: 
No Failure: 


u " u b * C P< Z P - W ♦ C I / < Z I - z ib> at 
Failure : 

u = u b + C P< Z P - z Pb> + C I / < Z I - z Ib> dt 


( 1 . 2 ) 


(1-3) 
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NO FAILURE 


Figure 1.6 OTA Algorithm Revision 2 [9] 
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The symbols are the same as in (l.l). Note that the estimates are being used 
in the integral trim portion after failure isolation. 

As part of the revision evaluation process, SCT: 

(1) generated new gain matrices for both the normal and failure 
mode filters used In the DIA algorithm; 

(2) defined a minimum complexity DIA algorithm; 

(3) examined estimator gain sensitivity to flight condition; and 

(4) developed the revisions for correcting the bias problem in the 
estimator outputs. 

As a result of this study, the advanced DIA technique (Revision 2) was 
shown to have improved performance over the technique chosen in the first 
program. However, model errors limited the performance of the estimators and 
made determination of detection and isolation thresholds difficult. The 
present study attempts to alleviate these problems. 

The scope of the present research effort is to develop an analytic 
understanding of the problems by applying the tools of robust multivariable 
control theory to accommodate modeling errors. 

1-3.2 The Fundamental Issue: Model Uncertainty 

The previous section has Identified model uncertainty as the main source 
of problems in sensor DIA algorithm design. Model (system, plant) uncertainty 
refers to the uncertainty in the errors between the nominal model and the 
actual system. There are two generic uncertainty representations : structured 
and unstructured [13]. The former refers to model parameters which are 
uncertain. The latter refers to unmodeled dynamics which are also uncertain. 
Reduced-order modeling techniques, linearization about operating points, 
neglecting nonlinearities, etc., all result in contributions to either 
structured or unstructured uncertainty. Figure 1.7 illustrates how 
uncertainties can appear in a control system which includes sensor failure DIA 
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Figure 1.7 Block Diagram of Control System with Sensor Failure DIA 
Logic, Showing Sources of Uncertainty 
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logic. The "blocks" represent fixed devices or processors. The "clouds" 
represent the system uncertainty and can be broadly grouped into two classes: 

(1) uncertain (external) inputs 
r - reference commands 

d - environmental disturbances 
b - biases or drifts in a failed sensor 

(2) uncertain (internal) dynamics, defined by the following 
transfer function matrices: 

a u (s) - plant model errors 

a^(s) - sensor failures 

a (s), a (s) - control law and estimator reconfigurations from 
C E DIA 

Since model uncertainty is the source of most difficulties, an approach 
is needed to deal with it in an effective manner. Robust DIA design provides 
a means to do this. 

1.4 METHOD OF APPROACH 

The method of approach developed in this report is to use robust control 
theory for threshold selector analysis and robust DIA filter design (as 
described In Section 1.2). This requires the isolation of uncertainties as 
shown In Figure 1.7. The approach establishes quantitative statements about 
the Interrelation among performance, robustness, and system uncertainty. 

Figure 1.8 Illustrates how all the uncertainties can be separately grouped for 
analysis. The external inputs, such as commands, disturbances, and failure 
biases, enter the system from the "outside." The dynamic uncertainties, such 
as model errors, sensor failures, and DIA reconfigurations, are "inside" the 
system and function as a feedback loop around the "interconnection" system. 

The interconnection system maps the external and internal uncertainties into 
the outputs, i.e. tracking error and filter residuals. 
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Fault-Tolerant System (with Separately Grouped Uncertainties) 


Up until recent years, there has been no unified control design approach 
for such a system as shown in Figure 1.8. Research efforts [13-15], which are 
elaborations on the Small Gain Theorem [16,17] have established the 
mathematical framework for such an approach. Basically, the dynamic uncertain- 
ties ( A|y|, A c> a £ ) all propagate in a specific way so as to cause a 
quantifiable uncertainty about the map from the inputs (r,d,b) into the 
outputs (e,y), provided bounds can be found for the dynamic uncertainties. 
These bounds are obtainable from simple input/output system tests. 

Furthermore, if the nominal system model is linear, the bounds and subsequent 
input/output errors are fully representable in the frequency domain. For 
example, the effect of all dynamic uncertainties on the tracking error and the 
filter residuals can be represented by simple graphs. For good tracking with 
no sensor failures, the error response and the filter residuals should be 
small over all frequencies. On the other hand, if a failure occurs — either 
a bias/drift b or structure change & s — then the filter residual 
frequency signature should be dramatically different from the normal 
(unfailed) mode. Otherwise, detection is not possible. 

The approach utilizes the following recent advances [13-20] in control 
theory: 

(1) Uncertainty Propagation [15] 

The dynamic uncertainties, as shown in Figure 1.8, can be bounded in 
a sector in the frequency domain which then propagates through the 
interconnection system so that the system input/output map is also in a sector 
in the frequency domain. These sectors determine quantitatively the 
performance/robustness trade. 

(2) Internal Model Principle [18] 

The principle states that asymptotic tracking and disturbance 
rejection In the presence of plant uncertainty can only be achieved if the 
controlled system contains a replica (internal model) of the commands and 
disturbance signal generators. For example, tracking constant commands 
requires Integrators in the loop. 
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These control design advances, together with methods for frequency-shaped 
LQG [26] allow for design of either the OIA logic or control/estimator 
directly in the frequency domain. Thus, if detection Is desired for a 
particular failure, the compensation required can be seen in the frequency 
domain exactly . Simultaneously, one can determine the effect of the new 
compensation on tracking performance. A similar procedure can be used for 
adaptive design. 

These ideas are illustrated In the flowchart of Figure 1.9, which shows 
how the tools generated from the method of approach can be used in a design 
process. Notice that the modeling errors combine with the sensor failures and 
the errors from OIA reconf igurations to form the system dynamic uncertainty. 

It Is possible to determine a bound on this uncertainty. The control system 
(estimator/control law/DIA logic) is then evaluated in the frequency domain by 
propagating the dynamic uncertainty as outlined above. The evaluation process 
yields quantitative results which suggest how to modify or robustify the 
design. This is an iterative process. 

Notice that there are three "levels of design. Constant gain OIA filters 
provide a "first cut" design. Optimal thresholds for these filters can be 
determined as well as minimum size of detectable failures (Level 1). Note 
that these computations are a function of bound on model error, noise, and 
class of failures. If a higher level of performance Is required, a robust OIA 
design can be carried out as shown in the flowchart (Level 2). Still better 
performance is possible by making the OIA filters adaptive (Level 3). 

1.5 OUTLINE OF REPORT 

This chapter has discussed the background and problem statement for this 
research effort. Current technology approaches were reviewed. The 
fundamental problem to be dealt with in OIA logic design, to achieve better 
performance, has been identified as model uncertainty. The proposed solution 
has been discussed at a high level. The rest of this report contains the 
technical details of the method of approach. 
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Chapter II discusses the effect of model uncertainty on stability and 
performance. Both stability and new performance robustness measures are 
developed. The performance robustness measures are with respect to tracking 
and disturbance rejection properties of the system. The effect of model 
uncertainty on asymptotic tracking capability is discussed. It is pointed out 
that for adequate asymptotic performance, an internal model is necessary. 

This creates certain so-called "structurally robust blocking zeros" to 
guarantee asymptotic performance. 

Chapter III discusses the effect of model uncertainty on sensor failure 
DIA. The fundamental concept of the threshold selecto r is introduced. This 
represents a new, innovative tool for the analysis and synthesis of OIA 
algorithms. The threshold selector concept is illustrated both for a scalar 
example and for a model of a multivariable turbofan jet engine. A closed-form 
solution has been obtained for the scalar case, and a computer-aided design 
(CAD) approach was developed to solve the multivariable case. 

Chapter IV discusses the design of robust filters for sensor failure 
DIA. The robustness of the filter is due to the presence of an internal model 
(as discussed in Chapter II) and due to the frequency-shaping of an IQG cost 
functional. The formulation and mathematical details of the frequency-shaped 
filters are discussed in this chapter. Results from an example using a 
dynamic simulation of a turbofan jet engine are presented. 

Chapter V presents the results of an evaluation study comparing the 
performance of the proposed robust DIA scheme to one of the schemes developed 
in Ref. 9. It is found that the performance of the two techniques Is similar 
for hard failures, but the proposed scheme shows considerable Improvement in 
the case of soft failures. 

Chapter VI provides some concluding remarks and directions for future 
research . 

Appendix A discusses model uncertainty and contains a procedure along 
with results for generation of a bound on model error for a jet engine. 
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1.6 REMARKS ON NOTATION 


Standard notation from control theory literature has been used whenever 
possible. For example, time domain quantities are usually denoted by lower 
case, whereas the same quantity In the frequency domain has been denoted by 
upper case. However, there are several exceptions. Lower case b(t) is used 
to denote a bias in the time domain and b(s) denotes the same quantity in 
the frequency, and similarly for v(t) and v(s). There are other 
exceptions that are generally clear from the context. 
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II. MODEL UNCERTAINTY AND ITS EFFECT ON STABILITY AND PERFORMANCE 


Model uncertainty is the main problem in DIA algorithm design as 
discussed In Chapter I. How model uncertainty manifests itself affects what 
one can do to achieve a high performance DIA algorithm. Hence, Its effects on 
stability and performance need to be studied before addressing its effects on 
detection and robust filter design (to be discussed in Chapters III and IV, 
respectively. This chapter discusses the sources of model uncertainty in 
dynamic systems. Both unstructured and parameter uncertainty are considered. 
Details on model uncertainty and a procedure for computation of its bound for 
a jet engine example are contained in Appendix A. This chapter concentrates 
on the effects of model uncertainty on stability and performance properties of 
the system. New measures are defined for performance robustness analysis 
which are similar to the well-known stability robustness measures [13]. The 
effect of model error on asymptotic tracking is explored. It Is shown that an 
Internal model is required in the DIA filter to eliminate the effect of the 
biases on the outputs of the system asymptotically. The inclusion of the 
Internal model results In creation of certain structurally robust blocking 
transmission zeros to guarantee robustness. 

2.1 MODEL UNCERTAINTY 

A very natural way to determine model uncertainty is to perform an 
experiment which compares the model with data from the actual system (plant). 
If there is no error between the model and the plant, then one has perfect 
knowledge of the plant. Normally, this Is not the situation the error is 
non-zero and represents how close the model Is to the plant. For example, 
consider the simple experiment depicted in Figure 2.1(a), where e is the 
error between the perturbation output iy^ = y - y Q from the actual 
plant and the corresponding output 5y^ of the linearized model; 5u 
is the perturbation input to both the plant and linearized model. P Q ( S ) 
a transfer function matrix description of the model. 
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By applying sinusoidal test Inputs, one can experimentally obtain a 
frequency-dependent bound as shown in Figure 2.1(b). Often, a good model will 
be accurate at a certain frequency range, small 6(w)*, and less accurate 
at other frequencies, large 5(w). Figure 2.1(b) Is characteristic of 
unmodeled high-frequency phenomena, incorrectly modeled parameters, as well as 
some types of neglected nonlinearities. For jet engine models, the key 
problem is "DC mismatch," i.e. non-zero error at low frequency. This is the 
predominant cause of setpoint "hang-off" and estimator bias. 

Sources of this type of uncertainty may be quite diverse. Slowly 
drifting parameters in an otherwise perfectly known linear time invariant 
(LTI) plant could yield the same uncertainty description as a plant with 
unknown nonl inearities approximated by an LTI model.** The key feature of 
this uncertainty is that, although it is bounded, one does not know the 
structure. Following Ref. 2, this is called the unstructured uncertainty of 
the model P (s). Essentially, the unstructured uncertainty Indicates the 
accuracy of the model in the neighborhood of the equilibrium. Equivalently, 
the plant input/output behavior in the neighborhood of the equilibrium can be 
described by the linearized uncertain model 

P (s) = (I + a(s))P ($) (2 ‘ 

m o 

where a( s) is the LTI unstructured output-multiplicative uncertainty 
operator with transfer function matrix a(s). All that is known about 
A(s) is that it is stable, causal, and bounded by 

o[A(jw)] < 4(u>), <*> > 0 (2 ’ 2 

where a(*) denotes the maximum (upper) singular value. 


* 6(<j) representing model error is net to be confused with perturbation 

quantities such as AYp, 6Y m , etc. 

** To reflect parameter errors in A(w) may require a very large number of 
experiments as shown in Figure 2.1(a). 
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In a general way, unstructured uncertainties account for all neglected 
dynamics, approximated nonlinearities, etc. 

One further remark about the uncertain linearized model P (s). For 

m 

every Input/output pair (u, y) which satisfies 

Y(s) = P(s) U(s) (2.3) 

in the neighborhood of the equilibrium, there exists a p m ( s ) 

(equivalently, a A(s)) such that Y(s) = P(s) U(s). One way to view this 
is to imagine that a(s) contains a sufficiently large (theoretically 
infinite) number of adjustable parameters so that any input/output can be 
perfectly matched. The bound on A(s) in (2.2) conveys the worst case 
situation. Note that a(s) can be infinite dimensional, without loss of 
generality. 

Parameter Uncertainty 


Modeling of plant uncertainty by the description given in (2.1) and (2.2) 
will be used throughout. A state-space linear plant model is assumed to be 
given by 


x = Ax + 8 u 
o o 


y = C x + 0 u 
oo 


with a corresponding model transfer function matrix 


P (s) = C (si - AJ 1 B„ + D„ 
0 0 0 0 0 


(2.4) 


(2.5) 


This model is valid only for Input/output behavior restricted to a 

neighborhood of the equilibrium point. A new equilibrium point will have a 

different nominal model, i.e., the parameters in P (s) will change as a 

o 

function of the equilibrium. This is easily seen in the definition of 
P q (s) given by (2.5). In a general way, one can consider the nominal model 
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to be a parametric model, where the parameters are adjusted to best fit the 
Input/output data. 

For example, let P (s) denote a parametric model of an uncertain 
plant P(s), where a is a k-vector of parameters in the model. Standard 
parameter identification methods can be used to find the best a c R to 
fit the data from the actual plant. For example, given plant input/output 
data (u, y) on the Interval t c[0, T'], a good parametric model Is found 

from 

inf ii y - y a ii T (2.6) 

acR k 

where inf denotes the greatest lower bound or infimum and IMI is a 
suitable norm, such as the norm 

II x ll T =^ T x 1 ( t ) x(t) dtj 1/2 < 2 - 7 ) 

The perfect matching condition 

Y (s) = P (s) U(s) (2,8) 

a a 

for some a and all (u, y) is never achieved for models of actual systems. 
The usual situation is the opposite. In fact, there is usually a range of 
a which solve (2.6) . 

The variations in a can be considered as uncertain parameters in the 
model. Since these parameters enter into the model in a definite manner, they 
can be referred to as the structured uncertainty of the model [13-15]. In a 
general way, unstructured uncertainties account for all neglected dynamics, 
approximated nonlinearities, etc. The structured uncertainty essentially 
yields the range of parameter variation in the model for a best fit of the 
data . 

Parameter (or structured) uncertainties can also be viewed in a slightly 
different way. Consider the system 
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(2.9) 


X ■= (A q + 4A) x + (B q + <$B) u 

y = C x 
o 

where (A q , B q , C q ) are the nominal system matrices and (6A, 4B) 
represents a perturbation in the system parameters. Suppose further that 
(4A, 4B) are known to be bounded such that 

o(6A) < a, o(4B) < b (2.10) 

Following the procedures developed in Refs. 4 and 5, one can construct bounds 
which include these parameter perturbations. Such a construction is useful 
when analyzing the effect of specific parameter uncertainties. Bounds on 
model error constructed in this way account for the type of uncertainties 
discussed by Leininger [4], 

2.2 EFFECT OF UNCERTAINTY ON STABILITY AND PERFORMANCE 

The model uncertainty discussed in the previous section has significant 
effects on both stability and performance of the system under feedback 
control. This section establishes quantitative trade-offs between uncertainty 
and stability, as well as performance. For this purpose, generic DIA 
configurations need to be considered for analysis and design purposes. In 
previous programs [8, 9], an estimator was designed as part of the DIA logic 
to provide synthesized estimates of system outputs. The estimator may operate 
both as part of the feedback loop or out of the feedback loop. 

Consider the two generic candidate design schemes shown in Figures 2.2 

and 2.3. The design in Figure 2.2, referred to as DIA 1, shows a feedback 

control system with an estimator running out of the control loop, i.e., 

"piggyback." The design In Figure 2.3, referred to as DIA 2, shows the 

estimator running In the loop. In both cases, the estimators contain the 

nominal plant model P (s), which approximates the actual plant P (s). 

o m 

and a dynamic filter gain F(s), such that 
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( 2 . 11 ) 


Y(s) = P Q (s)U(s) + F(s)v( s) 
v(t) = z(t) - y(t) (2.12) 

z(t) - b(t) + y ( t ) (2.13) 

where y(t) is the estimate of the output y(t), v(t) is the innovations 
signal, and z(t) Is the sensor output, which Is composed of the engine 
output y and a bias signal b(t) which represents both sensor noise and a 
class of sensor failures. The estimator structure (2.11) is equivalent to the 
usual state space structure for estimators/observers, where 


x = Ax +• Bu + Kv( t) , v(t) = y(t) - y(t) (2.14) 

y = Cx + Du (2.15) 

In this case, the nominal plant model and filter are, respectively: 

P Q (s) = C ( s I - A) _1 B + D (2.16) 

F(s) = C(SI-A) _1 K (2.17) 

The control signal is given by 

( G ( s) ( R( S)-Z(s) ) , DIA 1 (2.18) 

U(s) = 

( G (s) (R( S)-Y(S) ) DIA 2 (2.19) 


where R is the reference command and G (s) is the controller transfer 

c 

function. 

Both of these designs capture the significant characteristics of current 
DIA schemes as well as being general enough to account for a large class of 
alternate DIA schemes. For example, DIA 1 is representative of either: 

(1) a nominal (no failure) controller with DIA logic based on a 
"piggyback" estimator; or 

(2) a reconfigured system where G c ($) is a control ler/estimator 
combination with an estimator out of the loop for further OIA 
action . 
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Note that in (1) and (2) above, the out-of-loop estimator need not be the 
same. Similarly, DIA 2 is representative of either: 

(3) a nominal (no failure) system where the DIA logic Is based on 
an estimator in the loop; or 

(4) a reconfigured system where the DIA logic Is based on the 
estimator in the loop. 

As in (1) and (2), the estimators In (3) and (4) need not be identical. 

2.2.1 Effect of Model Error on Stability 

In this section, tools are developed for determining the effect of 
uncertainty on stability of DIA 1 and DIA 2. It will be assumed here that the 

plant is represented by 

P (s) * ( I + &(s)) P (s) (2 ’ 20 

m o 


with 


a[A(jw)j < i(w) , w > 0 

where P (s) Is the nominal engine model obtained by linearization about an 
equilibrium and A(s) is the unstructured model uncertainty. 


In order to analyze the systems, it is necessary to determine the 
transfer functions defined implicitly by, 

E ( s ) = H R(s) + H (s)b(s) 
er eb 

v(s) = H (s)R(s) +■ H (s)b(s) 

vP VD 

The most convenient way to do this is via the "interconnection" structure 
where the uncertainty A(s) is Isolated [15] as shown in Figure 2.4, 


E(s) 


H 

er 

H K 

eb 


’ R(s) 

v(s) 

J 


H 

L v> r 

H h 
vb j 


b(s) 


33 





A more detailed discussion of the interconnection system is provided in the 
next section. For DIA 1, one gets (suppress the 's' variable) 


H er = S( I + AT) 

H.=T+S(I+ AT) _1 AT 
eb 

H = L(I + AT) _1 AT 
\>P 

H . = L(I + AT)" 1 
vD 

and likewise for DIA 2, one gets 

H = S - S(I + GL) ( I + ATM) _1 AT 
er 

H . = TM + S(I + GL) ( I + ATM)’ 1 ATM 
eb 

H = L(I +• ATM)’ 1 AT 
vP 

H . = L(I + ATM)" 1 

vD 

where 

G = P G 
o c 

S = (I + G) -1 , T = (I + G) _1 G 
L = (I + F) -1 , M = (I + F) V 

The transfer functions (2.25) and (2.26) are stable if 

(1 ) S, T, M, and L are stable 

(2) 6( <j) <j[ T( ju) ] <1, w > 0 

(3) 4(«)o[T(jo»)M(ju)] < 1 , ^ > 0 

Consider the following scalar example, where 


(2.25) 


(2.26) 


(2.27) 


(2.28) 
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(2.29) 


P 

o 


(s) 


a 

s+a ’ 


F(S) 


k 

s+a 


G c (s) 



(2.30) 


and from (2.27) , 

r , i s = -L_ T = _a_ 
s’ s+a * s+a 

i = s+a M = k - 

s+a+k ’ s+a+k 


(2.31) 


Consequently, conditions (2) and (3) above yield the following model error 
bounds : 


( 2 1 ) «(«) < (l + (J) 2 ) 17 ? « > 0 (2.32) 

(3') «(«) < (l + (J) 2 ) 1/2 [0 + J ) 2 + (fc) 2 ] 1/2 . « > 0 


The above bounds can be interpreted as the maximum permissible bounds on 
model error for which the system (OIA 1 or DIA 2) remains stable. For 
example, if the actual model error is predominantly DC mismatch where the DC 
gain is known to within +10 percent, i.e., 5 = .1, then the above is 

certainly satisfied. 


2.2.2 Effect of Model Error on Performance: Performance Robustness 


Performance refers to the behavior of the system in relation to specified 
objectives, such as transient response, command following, and disturbance 
rejection. Performance robustness is the ability of the feedback system to 
maintain a performance specification despite plant uncertainty. In order to 
realize such robustness properties of feedback, it is necessary to establish 
quantitatively the trade-offs and relationships among performance, robustness, 
and plant uncertainty. Presented below is a method for directly analyzing 
performance robustness of uncertain multivariable systems. 
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2.2.3 Conic Sectors 


Structured and unstructured uncertainties can both be viewed as belonging 

to an L -conic sector, defined as follows [15]: 

P 

Definition : A relation (u, y) c H is inside the Laconic sector, 

denoted H c L -Cone (C, R, S) if for some p c [1 , «] there exists 

compatible ope?ators C, R, and S such that H - C is testable* and 

nS(y - Cu) 11^ < nRuiip for all (u, y) c L p x L p with y c Hu. 

In the case of Instable LTI operators R, S, and H-C with transfer 
function matrices R(s), S(s), and H(s)-C(s), respectively, the conic 


sector is equivalent to the frequency-domain condition. 


tf[S(j«)(H(j«)-C(J»))R 1 (j«)] 1 1. 
Let n = diag (a 


u 


(2.33) 


a ) be a diagonal matrix containing all the 
k 


structured uncertainties. Equivalent statements for plant uncertainty with the LTI 
model P [a] of P is that the structured uncertainty Q c L^cone (0, 13, I), 
where (3™= diag (13 ... 8.) and the unstructured uncertainty r c L^cone (0, 

I, I). Thus, conic sectors conveniently describe typical model (plant) uncertainty. 


2.2.4 Plant Interconnection Model 


It is convenient to represent the uncertain plant as shown in Figure 2.5, 
which Is patterned after the uncertain system descriptions presented in 
Ref. 15. The operator 
Q 0 


(2.34) 


* A relation H is said to be L p -stable if for all u e L p , y e Hu c L 
and there exists finite positive constant k such that llyii p < kilull p . 
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contains all the uncertainties in the plant — both structured and 
unstructured. The operator H Is referred to as the plant interconnection 
system and serves to Isolate all the uncertainties in a from the rest of 
the plant. Consequently, H is completely known. The governing plant 
equations are then 



V = AZ 

The plant input/output relation Is given by 
y = P[A]u 


with 

P[A] = P + H (I - AH ) _1 AH 

o yv zv zu 

where P Is the nominal plant defined by 
o 


(2.35) 

(2.36) 

(2.37) 

(2.38) 

(2.39) 


The uncertainty A can also be viewed as belonging to a conic sector, 
i.e. a c L^-Cone (0, R, I), R = d1ag(0, I) which implies the conic 
sector descriptions of the uncertainties of Q and r. The plant 
representation shown In Figure 2.5 is extremely useful in evaluating 
performance of the closed-loop system. With an LTI model, the conic-sector 
bound on A is equivalent to the frequency-domain condition. 


a(A(jc J )R _1 (j w )] < 1, 


R(ju) = diag(l3, l( w )I) 


(2.40) 


2.2.5 Normal i zation 

The uncertainty and interconnection system can be normalized . Let 
AR ^ replace A, RH zv replace H zv> and RH zu replace 


38 




39 




P[fl] now has the same form as before, but a Is normalized so 


H zu’ 

that 


a c Lp-cone ( 0 , I, I) 

For convenience, assume that henceforth these replacements have been made and 
the plant interconnection system is normalized. 

2.2.6 Performance Robustness Measures 


Consider the feedback system shown In Figure 2.6. The plant has the 
transfer function 


P(&) « P Q (s) (I +- a(s) ) 


(2.41) 


where P (s) is the nominal plant transfer function and a(s) represents 
o 

the unstructured input-multiplicative uncertainty. P (s) expressed in 

o 

terms of state variable matrices is 


P 0 (s) = C(sl-A)" 1 B f D 


The closed-loop interconnection system is 
Y 


YR (S > 

H Vd (S) 

h yv ( s > 

ZR (S) 

H Zd< !) 

( S ) 


R 
d 
L V 


where for a unity feedback system ("suppressing" s) 


-1 


YR 


( I+P o G c ) P o G c 


H YV=< I+P o G c> P o 


’ZV 


(I 4- G P ) 1 G P n 

CO CO 


-1 


H ZR = ( I+G c P o> G c 


(2.42) 


(2.43) 


(2.44) 

(2.45) 

(2.46) 

(2.47) 
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UNCERTAIN PLANT 
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Figure 2.6 Control System with Uncertain Plant Feedback 





Now assume that the tracking performance requirements are expressed in the 
frequency domain in terms of 


"V a) - y 

HH yR ll 


< p(u) 


co > 0 


(2.48) 


where H v _ Is the nominal input-output transfer function, H ( A) is 
the perturbed input-output transfer function and p(«) represents some 
given function of frequency which specifies tolerable tracking performance 
degradation. Furthermore, if 4(w) represents a bound on model 
uncertainty, i.e. 


HAD < 6(u) co > 0 


(2.49) 


then the tracking performance robustness properties of the system Is described 
by the following theorem. 


Theorem 2.1 : Assume that the tracking performance requirements are 

expressed by (2.48) and the closed-loop interconnection system is stable. 

Then tracking performance requirements are guaranteed if 

J (u) > <S(u) CO > 0 (2.50) 

pr t 

where 


p(u) O(Hyp) 

5 pp (w) = (2.51) 

T a(H zv )c(H YR )p(<o) 4- 0 (H yv )o(H zr ) 

is the tracking performance robustness measure. 

Proof : Assuming d = 0, from Eqs. (2.38) and (2.43), it follows that 

H yr (A) = H yR +■ H yv (I - AH ZV ) ] A H zr (2.52) 

Substituting (2.52) into (2.48), one obtains 

HH YV ( I “ AH ZV )_1 A H Zr" - ,iH yr" (2<53) 
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Using the matrix inequalities 


o(AB) > cr( A) o(B) 

o(I - C) > 1 - o(C) if o(C) < 1 

3(A) = — -T- 
o( A ') 

then it follows from (2.53) that 


(2.54) 

(2.55) 

(2.56) 


„<H yv ) 1 o(H 2R ) _ 

• < P(w) 

1 - 5 o(H zv ) 


«(h yr ) 


(2.57) 


Solving for the worst-case value of 6 in (2.57) and naming the particular 
value satisfying (2.57) 4 pR we have 


5 


PR, 


p(w) '’(Hyp) 


(2.58) 


Therefore, so long as 

4 < 4 pr t ‘ 

then (2.57) is satisfied and the tracking performance robustness is guaran 
teed. Q.E.D. 

Suppose that the disturbance rejection performance robustness is 
expressed as 

HYll < 0(o,) w > 0 ( 


where 0(u) describes the allowable disturbance rejection performance 
degradation. The following theorem describes the performance robustness 
properties of the system. 
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Theorem 2.2: Assume that the disturbance rejection requirements are 

expressed by (2.60) and the closed-loop interconnection system is stable. 

Then the disturbance rejection performance robustness is guaranteed if 

( w ) > £(<*>) a) > 0 (2.61) 

where 



is the disturbance 


o(H zv )(B-o(H Yd )) + 0(H yv MH Zd ) 
rejection performance measure. 


(2.62) 


Proof: Assuming R = 0, from (2.36) and (2.43) it follows that 

H Yd (A) = H Yd * H YV fl(I " H ZV A)_1 H Zd (2.63) 

and 

llH Yd + V^ 1 ~ H ZV A ^ H Zd M < B (2.64) 

Using the matrix inequalities (2.54) through (2.56), the above equation can be 
written as 

«(H Yd )(l - o(H zv ) 6) + 0(H) 5 o(H ) 

I < B (2.65) 

1 - 6 0(H ZV ) 

Solving for the worst-case value of <S from (2.65) and naming it 6 pR («) 
we obtain ® 


, , , B - 

^PR ~ I “ 

D 0(H zv )(B - 0(H yd )) + 0(H yv )5(H Zd ) 


( 2 . 66 ) 


Then, so long as 
6 < 6 


PR 


D 


(2.67) 


(2.65) is satisfied and the disturbance rejection performance robustness is 
guaranteed. Q.E.D. 
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( 2 . 68 ) 


The stability robustness measure for the system is given by [13]: 


4 SR (<j) 


1 

^(H zv ) 


and the system is stable provided that 


(2.69) 

(2.70) 

(2.71) 


6 sr (“) > *(«) 

It is easy to see that 

‘pr/” 1 ' 4 SR M " - 0 

‘rr^ 1 < J SR ( “’ “ 1 ° 

Figure 2.7 shows the stability and performance robustness measures for a 
fifth-order multivariable system for 5% performance degradation from nominal. 
These measures then establish how accurate the plant model has to be for 
specified performance as a function of frequency. Note that various 
controller designs can be evaluated on the basis of these measures. 


2.2.7 Effect of Model Error on Asymptotic Tracking 

The discussion to follow will Illustrate the Internal model principle 
[18], namely that, design of robust filters (and controllers) to achieve 
asymptotic tracking and disturbance rejection, despite model error, can only 
be accomplished if a model of the command and disturbance Is incorporated In 
the filter (controller). It is shown that the Inclusion of the Internal model 
creates appropriate structural ly robust blocking zeros. 

The critical performance measure for engine control is the ability to 
achieve asymptotic tracking despite model error at OC. This section 
highlights some of the main issues associated with this problem. 

In order to evaluate the effect of model error on asymptotic tracking, 
assume that the following conditions hold: 
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(1) The model error satisfies (2.28), i.e., the system Is stable. 

(2) The reference r(t) -» r ss (constant) as t -* ®. 

(3) There is no sensor failure, i.e., b = 0. 

These conditions, together with (2.31), yield the following steady-state 
signals for OIA 1 (estimator out of loop): 

e = 0 (2.72) 

ss 


H vr = L( I+AT) ] AT 


s+a . _a_. -1 . _a_ 

s+a+k + a s+a s+a 

(2.73) 

(s) = H (s) R(s) 
\>r 

(2.74) 

m _ MIO) R (s) 

ss iSJ ' (a+k) ( 1 +A(0) ) R ss Vb; 

(2.75) 


Similarly, for OIA 2 (estimator in the loop) 


v(S) = H yr (s)R(s) 

(2.76) 

H = L(I+4TM) _1 4T 

vT 

(2.77) 

aa(0) r 

e ss v ss a^-k^kA(O) ss 

(2.78) 


This latter result exhibits the undesirable "hang-off" when the estimator is 
In the loop, e.g., OIA 2. This problem can be eliminated if the filter in 
OIA 2 contains an integrator, e.g., if F(s) In (2.30) has the form. 


F(s) 


f 1 s +f 2 
s(s+f 3 ) 


(2.79) 


To see this, consider the plant and controller of (2.29) with F(s) as in 
(2.79), then 


G 


a 

s 



(2.80) 
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L = (I+F) 


-1 




s(s + f 3 ) + f ]S + f 2 


(2.81) 


M = IF = 


V * [2 


5(S+f 3 ) + f ^ S + 


(2.82) 


For OIA 2 


H = L(I + ATM) AT 
vf 


(2.83) 


and 


11m H - 0 

a vr 

s-*0 


(2.84) 


It then follows that e = 0. 

ss 


Internal Model Principle and Structurally Robust Blocking Zeros : 


The internal model principle states that asymptotic tracking and 
disturbance rejection in the presence of plant uncertainty can only be 
achieved if the controlled system contains a replica (internal model) of the 
commands and disturbance signal generators [18-20]. For example, tracking 
constant commands and rejection of constant disturbances requires integrators 
in the loop. The presence of the internal model creates zeros of transmission 
In certain transfer functions. For instance, for command following, transmis- 
sion zeros are created between the reference input and the tracking error at 
the frequencies of the commanded signal. These zeros referred to as 
"blocking" zeros, are also structurally robust , i.e., the locations of the 
zeros are unaffected by parameter variations. 


For sensor failure DIA systems, if the plant's steady-state output is to 
be unaffected by a failure (represented by the bias b), then the transfer 
function from b to u (and therefore y) must contain a structurally robust 
blocking zero. This would happen if an internal model is present in the 
filter. The following examples illustrates this. 
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Example 2.1 SISO System with Unknown Measurement Bias: 


Consider the SISO plant 

p 0 (s) = 7 ^ 

with the filter containing the internal model 
Plant: 

x = -ax + au 
y = x 
z = y + b 

■ A 

Filter: x = -ax + K^z - z) + au 


(2.85) 


( 2 . 86 ) 


b 


K 2 (z " z) 


(2.87) 


z = x + b 
Control law: 

u = -K q x 


( 2 . 88 ) 


Then substitution of the control law in the system equations results In the 
overall system equations from the bias to the control 






- 


- - 



• 






- 



X 


-a- K r aK o 

- K 1 

K 1 


X 


K 1 

ft 






- 



b 

= 

- K 2 

' K 2 

K 2 


b 

+ 

! 

K 2 

• 

X _ 


--* K n 

0 

-a J 


- x J 


.0 . 


u = -[K q 0 



x 


(2.89) 
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For a linear time-invariant system of the form 


x = Ax + Bu 
y = Cx + Du 


( 2 . 90 ) 


the transmission zeros are defined as those frequencies at which the system 
matrix (pencil) 


Y z (s) - 


S I -A 
-C 


loses rank [22]. 

The zeros from b to u are given by 




s+a+K. +aK 
1 0 

K 1 

- K 1 

1 

1 

K 1 

Y Z (S) = 

det 

K 2 

S+K 2 

- K 2 

1 

1 

K 2 



aK 

0 

0 

s+a 

1 

1 

0 



1 

O 

0 

0 

1 

1 

0 


( 2 . 91 ) 


= -K K.s(sfa) 
o 1 


( 2 . 92 ) 


i.e., there is a zero at the origin. Furthermore, this zero Is structural 1 v 
robust , i.e., its location does not change in spite of the parameter errors in 
the system matrices. To see this, assume that the system matrices are modified 
such that 


P 0 (s) - 


a+P 

s-t-a+c 


( 2 . 93 ) 


where c and P represent plant parameter variations of any size so long 
as the closed-loop system remains stable. 


The overall system equations from b to u are now given by 
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ft 

X 


(-a-aK 0 -K 1 ) 

- K 1 

K 1 


X 


K i 

X 

b 

• 

X 

= 

- K 2 

-(a+B)K Q 

- K 2 

0 

K 2 

-(a+c) 


b 

X 

+ 

*2 

0 

— 


X 


U = [-K, 


0] 


(2.94) 


The zeros are 


Y z (s) = det 


12 ' 


— 


- K ! 

1 


s+a+aK +K, 
o 1 

s 

1 

K 1 

K 2 

S+K 2 

-K 2 

1 

k 2 

(a+B)K 

0 

0 

s+a+e 

1 

1 

0 

K 

0 

a-c) 

0 

0 

1 

1 

0 


(2.95) 


wnicn confirms the fact that the zero at s=0 is unaffected by parameter 
variations hence the robustness property. 


Fxamnle 2.2 Multivariable System with Unknown Constant Measurement Bias 


Consider the general square multivariable system 

x = Ax + Bu 
Plant: y = Cx 

z = y + b 

* « * * 

Filter: x = Ax + Bu +• K^z-z), y = Cx 

z = y + b 

b = K 2 (z-z) 


(2.96) 


(2.97) 


(2.98) 
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Control law: u = -K x 

o 

The overall system equations are 


X 






— 




X 


A-BK -K.C 
o 1 

" K 1 

K-jC 


X 


K 1 

* 









b 

— 

-k 2 c 

- K 2 

K ? C 


b 

+ 

CM 

v: 

• 

X 


_ " 8K ° 

0 

A 


X 


0 


u = [ -K 0 
o 


0 ] 


transmission 

zeros of the 

system are 

given 

by 



Sl-A+BK tLC 
o 1 

K 1 

-K, c 

1 

1 


Y z (s) = det 

V 

si+k 2 

-k 2 c 

1 

1 



-B K 

0 

0 

Sl-A 

l 

1 

0 


-K 

o 

0 

0 

1 

0 





1 

— 


= 0 


which after elementary row and column operations is equivalent to 


y z (s) 



si -A 

" K 1 

0 

i 

1 

1 

I 

0~ 

det 

-K 

0 

1 

1 

1 

0 

1 

0 


0 

0 

1 

si -A 1 

0 


0 

~ K 2 

1 

1 

“I ‘ 

o 1 

Si 


which shows the presence of blocking zeros at the origin. To see that 
transmission zeros are structurally robust, consider the parameter and 
controller perturbations 

A = A + 4A 


B = 8 -i- 6B 


( 2 . 99 ) 


( 2 . 100 ) 


( 2 . 101 ) 


( 2 . 102 ) 

these 

( 2 . 103 ) 

( 2 . 104 ) 
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(2.105) 


K = K + 6K 
0 0 0 


K 1= K-j + 6K ] 


K 2 = K 2 + 6K 2 


The transmission zeros of the perturbed system are given by 


y z (s) = 


si -A 

“*1 

1 1 
1 o | 

0 

~ _ 

-K 

0 

0 

0 

o 1 l 

0 

0 

1 1 
</> 1 

0 

0 

-*2 

1 1 
1 ° 1 

si 




1 i 

— 


(2.106) 

(2.107) 


(2.108) 


which confirms the fact that the transmission zeros at the origin are 
structurally robust as they are not affected by plant parameter variations 


The above development can also be carried out for systems with a non-zero 
direct transmission term (D * 0) and with more complicated dynamics for the 
bias b. The conclusions would be the same. 

2.3 SUMMARY 

In this chapter, we have discussed various sources of model uncertainty 
and its representation. The effects of model uncertainty on stability and 
performance was discussed for two generic configurations. New performance 
robustness measures were Introduced to determine tracking and disturbance 
rejection performance robustness properties. The effects of modeling errors 
on asymptotic tracking were discussed. It was shown that an internal model is 
required in the filter to produce unbiased estimates. 
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III. EFFECT OF MODEL UNCERTAINTY ON FAILURE DETECTION: 
THE THRESHOLD SELECTOR 


One of the most difficult problems in sensor failure OIA algorithm design 
is the ability to evaluate analytically the DIA algorithm in the presence of 
model error. In this section, for the first time, a unified framework which 
allows for this analysis is presented. A new concept is introduced, referred 
to as the threshold selector , which is a nonlinear Inequality whose solution 
defines the set of detectable sensor failure signals. The threshold selector 
is consistent with the frequency domain model uncertainty description that has 
been emphasized in this study. What follows is a heuristic discussion of 
failure detection which leads to the notion of the threshold selector. For 
illustrative purposes, the focus will be on the innovations approach to 
failure detection. As will be seen, the methodology is quite general and not 
limited to just the innovations approach. The conventional technique of 
selecting thresholds in innovations-based DIA filters has been based on 
noise. In a previous study [9], constant thresholds were selected based on a 
measure or innovations size obtained from a no failure hypothesis filter. 
However, the current technique determines thresholds based on model error, 
sensor noise, and class of failures, as well as the speed of the DIA filters. 

The threshold selector inequality to be presented here represents a new 
and innovative tool in the analysis and synthesis of DIA algorithms. In 
particular, good estimates for the minimum threshold set are obtained in the 
multivariable case. In this case, it is necessary to compute operator gains 
dependent upon the norm measure used in detection. DIA designs, more sophis- 
ticated than those illustrated here, can also be analyzed and synthesized 
using the threshold selector inequality. 

3.1 INNOVATIONS APPROACH TO FAILURE DETECTION 

Many sensor failure detection schemes (for example, see Refs. 1, 8 and 9) 
have been based on determining the characteristics, such as the RMS value, of 
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the innovations v(t), over a given finite time interval t e[0, T]. In 
this section, a framework will be presented within which it is possible to 
evaluate directly the effect of model uncertainty on the ability to detect a 
failure. The analysis presented also proves to be extremely valuable in the 
selection of the 0 1 A estimator. 

The problem formulation is as In Chapter II. From (2.24), the Innovation 
sequence may be expressed as 

v(s) = H (s)b(s) + H (s)R(s) (3J) 

vD vr 

where H and H are dependent on the model uncertainty A, and from 

vb vr 

(2.25) and (2.26), 

H (s) = H (s)A(s)T(s) (3 ’ 2) 

vr vb 

Define the bias b as 


where n is zero-mean sensor noise and f is a drift signal associated vith 3 
class of sensor failures. In practice, the situation is that a known class of 
failures is possible, and detection is limited by some bounded noise signal 
and bounded model error. Let Q^, , and denote, respectively, 

the bounded sets of noise signals n, failure signals f, and model errors 
a. Let J( T ) denote a measure of the innovation size on the Interval 
0 < t < t, where J(t) is given by 

J( T ) = HvIIt 

The (truncated) norm operation imi is based on an RMS measure of the 
innovations and will be defined precisely below. The norm in Eq. (3.4) may be 
evaluated either in the time domain or frequency domain using Parseval's 
theorem [28]. Hence, the notation 

J(t) = Hv( s) II (3 ' 5) 
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denotes evaluation In the frequency domain. Substituting (3.2) into (3.4) 
gives 

J(t) = II H (s)(N(s) + F ( s ) + W ( s ) T ( s ) R ( S ) ) II (3.6) 

VD X 

With no model error, a = 0 and 

v($) = H yb (s)(N(s) + F( s ) ) (3.7) 

Thus, a non-zero bias due to sensor failure is detectable in y through the 
dynamics of H^ b given by (2.26). Under mild constraints, it is possible 
to detect the presence of failure for relatively small bias. The only limit 
is due to noise levels in the sensor. However, the model error effectively 
raises the detection threshold (3.6). It was found in previous applications 
of sensor failure DIA logic to the turbofan engine problem [8, 9] that the 
model error considerably dominates any sensor noise level and is the primary 
cause of false alarms and misses. What is sought is a means to detect the 
presence of failure bias b despite model uncertainty a. More 
specifically, what are the conditions of the transfer function matrices 
^yb an< ^ suc ^ that the ability to detect a failure is 

maximized? Further, what is the most sensitive detection scheme for a known 
class of failures? 

3.1.1 False Alarm 


A major requirement on detection is to reduce or prevent false alarms. 
Thus, In the absence of any failure signal, J(-r) should be less than a 
threshold value, J . Setting f = 0 in (3.6) gives* 

J th (r) = sup llH vb (s)(N(s) + W(s)T(s)R(s))ll (3.8) 


* 


For the sake of brevity, the notation sup(«) is used in place of 

x 

for all x = n, f, a. Likewise, sup(*) means sup(*) 

x.y xefl x ,ycfl x 


sup( •) 
X£fl x 
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where sup(*) denote the supremum or the least upper bound. Note that the 
threshold function Is dependent upon the known reference command r as well 
as bounds on the noise and model errors. In fact, with no model errors, 
a * 0 and the threshold is determined strictly by the worst-case noise 
level, i.e., 

J th (x) = sup llH vb (s)N(s)ll T (3.9) 


3.1.2 The Threshold Selector 


Setting a threshold to eliminate false alarms due to sensor noise and 
model error will cause the detector to miss certain failures. The question 
Is: What is the minimum detectable failure? In other words, find the 

threshold failure set. Denote this set by which is defined as 

those fcQj. such that 

inf J( T ) > J Jt) (3.10) 

. c th 

A,n,f 

where i nf ( * ) denotes the infinum or the greatest lower bound. Substitut- 
ing (3.6) and (3.8) into the above gives 


inf IIH (s)(N(s) + F(s) + a(s)T(s)R(s)il > sup iiH. (s)(N(s) + a(s)T(s)R(s)n 
a.n.f ub T a,n vb 

(3.11) 

If we decompose the Innovations Into components due to the noise, failure, and 
error model 


v(s) = v (s) + v (s) + v (s) 
n f a 


then (3.11) may be rewritten as 


inf iiv + v, + v II > sup iiv + v II 

AnC n f A T . „ n A T 

a, n , f a, n 


(3.12) 


(3.13) 


Since this inequality generates the minimum threshold set, it is referred to 
as the threshold selector . An estimate of the smallest size of failure f 
which is detectable can now be calculated. 
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Theorem 3.1 : An estimate of the size of minimum detectable failure 

| f | Is given by 

Ifl > 2 J th (t)/B(t) (3.14) 

with the threshold as 


J th (x) = (max o [ L( ju>) ] ) n + 6|r|o L(Q t1 LT)(t)] (3.15) 


and 


3(t) = o[Q t1 L(t - t f ) ] (3.16) 

| r | = norm of the reference input signal 
|f| = norm of the failure signal 

Q . = operator gains which are functions of type of failure and 
T reference signals; k = i,j 

6 = upper bound on model error (constant) 
t^ = failure time 

t = size of detection window 


n = bound on RMS sensor noise 

Proof : A conservative estimate of the threshold set can be found from 

the Inequality: 

Inf ll v, ( s ) ll > 2 sup II v (s) + v ( s ) ll (3.17) 

c*' T . _ n A t 

r , A A , n 

since (3.17) implies (3.13). This can be seen as follows. We may rewrite 
(3.17) as 
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jnf ii v f ( s ) ii T - su^ Hv n (s) + v A (s)ii > sug «v n (s) + v a^ s)ll T 

and since 

sup iiv (s) f v (s) ii > Inf Hv n (s) + v & (s)» T 
a.n n f - A 

( 3 . 18 ) may be written as 

inf Hv f ( s) H - inf Hv (s) * v (s)ii > sup iiv (s) * ^ a (s)i« t 
f , a f T a.n n 

However, 

II v f ( S ) + v n (S) + v 4 (S)ll t > |llv f (s)ll t - llv n (s) + v a (s)ll T l 
= llv f (S)ll T - llv n ( S) + V & (S)II T 


(3.18) 

(3.19) 


(3.20) 


(3.21) 


as long as 

|| v f (s)i | T > llv n (s) * v a (s)ii T (3l22) 

Suostituting (3.21) into (3.20), we have 

inf II v (s) + v f (s) + V (S)ll > sup II V (S) + v a (S)ll T 
a,n, f n f A - n 

(3.23) 


which is (3.13) . 

In order to utilize (3.13) or (3.14), it is necessary to specify the 
detection measure imi^, and the sets £2^., and Suppose that 

detection is based on the root mean square (RMS) measure, 


11X11 = C- f |x(t)| 2 dt) 1/2 

T T *'0 

where |x(t)| is the Euclidean norm of x(t), i.e 


1 /2 

| X ( t) I ■ (X' (t)x(t)) 


(3.24) 


(3.25) 
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Furthermore, let 


= {f I f ( t) = f Q l(t- t f ), | f Q | > f, t f c[0, r]}* 


(3.26) 


n = {n| iimiT < n ) 
n - 


(3.27) 


- (4|o[a(jw)] < 4(w), W > 0) (3.28) 

In other words, the sensor failure signal f occurs at some time t^e[0,x]; 

and Is represented by an abrupt shift in bias; the noise signal n is 

arbitrary, except that it Is bounded in norm (3.27); and the model error a 

is bounded as described in (3.28). Note that by allowing the noise to be In 

the set n , it is not possible to distinguish at each Instant of time n 
n 

cSi from f eSl . Therefore, it is necessary to view the innovations 
n f 

over a (moving) time window t, e.g., as in (3.24). One can now calculate 
an estimate of the smallest f in (3.26) by evaluating (3.13) over (3.17). 

The mathematical machinery for this calculation is available [14-16], but 
requires the introduction of the following definitions. 


Let H denote an operator with proper rational transfer function matrix 


H(s). Let 


denote the linear operator defined by 


<r k H)(t) = 


a *’ 1 
!£ - 1 


[H( s) ] , 
(k-D! 


H(s) 


k= 0 

k = 1,2,. 


(3.29) 


where '£ [•] is the Inverse Laplace transform operator. Thus, 

corresponding to H(s), (f Q H)(t) is the impulse response matrix, 

( r i H ) ( t ) is the step response matrix, and so on. It is also convenient 
to define the matrix operator. 


* l(t) is the unit step function, i.e., l(t) = 0, t < 0; l(t) =1, t > 0. 
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(3.30) 


(0 Tk H)(T) - [(r k HHt)]'[(r k H)(t)]dtj ’ /2 


With definitions (3.29) and (3.30), one can now calculate (3.17). In 
order to facilitate the presentation here, we make the simplifying assumption 
that the model error a is constant (e.g.. DC model mismatch only) and is 
sufficiently small so that 

(I + aT) -1 -I. (I + ATM) -1 ~ I (3- 3 


or equivalently* 

H h (s) = l(s) 


(3.32) 


From (3.17) it then follows that 

1nfnl(s)F(s)n > 2 sup nL(s)N(s) + L(s)a(s)T(s)R(s)M t 
f,a T a, n 


We will now evaluate the above norms In the time domain. 
[28] 

'laAll = |a| >1 A II 
II AB ll < It All II B II 
and the fact that 


Using the relations 


(3.34) 

(3.35) 


o(A) 


ll Ax ll 
- nxll 


< o( A) 


(3.36) 


the left-hand side of (3.33) may be replaced by 

I f I 2 [0 ,(L(t - t ))] (3 * 

t I r 

where the 1 n f ( * ) has been replaced by the lower singular value from (3.36) 
and (3.26) was substituted for f (* t ) . We then have that 


* Note that without this simplification, one has to deal with (3.17) directly 
which is rather awkward. 
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(3.38) 


I f I 2[0 ,(L(T-t ))] > 2 sup 11 L ( s)N(s) + L( s)W( s)T( s) R( s) II 

fl,n T 


Now If we ensure that 


|f| 2C0 .(Kt-t.)] > 2 sup( HL(s) II HN( s) II h- if A( s ) II nL(s)T(s)R(s) II 

T '' T T T T 

a , n 

(3.39) 

then using (3.35) and the fact that 

llA+Bll < II A II + ll B II (3.40) 

(3.37) follows from (3.38) since the right-hand side of (3.38) is smaller than 

(3.38) . Furthermore, if we use the identity [28] 


nail 2 = o(a) = max |A(jw)| (3.41) 

(O 

then the right-hand side of (3.39) may be written using (3.27), (3.41), 

(3.28), and (3.30) as 


2(max | l( jw) | n + 4|r|5[Q Tl (LT) (t) ]) 

Substituting (3.42) for the right-hand side of (3.39) we now have 
|fl2[Q T L( T -t f )] > 2 (max |L(j«)| n + 6 1 r | «CQ xl ( LT) < t ) ] ) 

1 u 


where we have used the assumption that 
r(t) = r l(t) 

If we define the threshold as 

J th (T) = max | L( j«) | n + 5|r| o[ 0 t1 ( LT) (t)] 

u 

and 

B(t) = 0 Tl (L(x-t f )) 
then (3.43) may be written as 


(3.42) 


(3.43) 


(3.44) 


(3.45) 


(3.46) 
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which is the desired relation (3.14). Q.E.D. 


Figure 3.1 shows a conceptual plot of |f| vs. t. Note that in the 
presence of noise, the detection window must be large enough to separate noise 
from the bias shift due to\ensor failure. Further, there Is a detection 
window t, dependent on t^, such that |f| Is a minimum. But, 
t is not known; consequently, one can only evaluate the effect of window 
selection. Figure 3.1 is also useful in evaluating the filter dynamics. If 
the dynamics of the filter are fast, there Is a sharper, lower threshold than 
when the filter Is slow. Figure 3.1 illustrates this critical trade-off In 
filter design, i.e., threshold vs. detection window. Figure 3.2 shows a 
typical plot of threshold J th vs. t. In Theorem 3.1, we assumed that 
the reference input signal was a step. However, other reference input signals 
may be used such as a ramp input signal, as done In the examples below. 

Example 3.1 : Consider the following scalar 


with proportional plus integral control 


example, where 


(3.48) 


g c (s) - 1*7 


and filter dynamics. 


F(s) 


s+a 


then from (2.27), 

G(s) ■ S(s) = s+g , T(s) = s+a 

= s+a+k ’ = s+a+k 


(3.49) 


(3.50) 


(3.51) 

(3.52) 



£( t): smallest size of detectable failure 
t, : time failure occurs 

t: length of detection window 


Figure 3.1 Threshold Selector: Minimum Detectable Bias Shift (f) vs. 

Detection Window (r) 
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THRESHOLD J 



J: Threshold 

J • Threshold for large detection window 
n: RMS noise 

t: length of moving detection window 


Figure 3.2 Threshold vs. Detection Window 
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and 


L < S > T < S > = I7ITk 


(3.53) 


r, LT(t) = 2*' 


-i a 


s ( s+a-t-k ) 


-(a+k)t 


r ] L(t) 


,y.-l s+a 

s(s+a+k) 


a k ~(a+k)t 
a+k + a+k e 


Therefore, 


Q ,i LT< ’> -(;/ 0 [* (1 - e ’ < atk ) t >] 2 « tV /2 




and after the Integrations have been performed, we obtain 


Q t1 LT(t) 


— ^ [2(a+k)T - 3 - e- 2(atk >' ♦ 4*-< atk b >' 2 

2(a+k) t 


Q Tl L (T-tf) 


r [2a 2 (a+k)(x-t f ) - k 2e- 2(a+k)(T - t f ) 

2x(a+k) 

, , -(a+k) ( -r-t ) ,2 . . . 1 1/2 

- 4ake f + k + 4ak] i 


(3.59) 
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Since the singular value of a scalar quantity is simply its absolute 
then , 

o[(Q t 1 LT)(t)] = |(Q xl LT)(t)| 

and 

o[(Q t] L) (x-t f )] = |(Q t1 L) (x-t f ) 1 
Also we have that 


<j(Uju)) = 


i 


? . 2 
3 (i) 


2 2 
(a+k) +• u 


and 


max o[L(jw)] = 1 since a > 0 and k > 0 

Finally, we obtain for the threshold 

J th ( T ) = n + 6 | r | |(Q t -| LT)(t) | 
or 


t ) = n + 5|r| • 



a 2 „ -2(a+k)x . -(a+k)T] 

— - — — [2(a+k)T - 3 - e £V ; + 4e 

1/2 


2(a+k) t J 



0(T) = 1(0 , L) (t) | 

T I 


value, 

(3.60) 

(3.61) 

(3.62) 

(3.63) 

(3.64) 

(3.65) 

( 3 . 66 ) 
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or 


B( 


T> = I ! 


2x(a+k) 


i - ? -2(a+k) (x-t ) 

J 3 [2a <: (a+k)(T-t f ) - ITe 


- 4aKe 


-(a+k)(r-t f ) 


+ K 2 + 4aK]i 1/2 


(3.67) 


The closed-form solution of threshold as a function of the detection window Is 
given by 


f(t) 


LV T) 

8(0 


In the limit as t + ® 


(3.68) 


J th (“) ■» " + 4 1 r | 


a 

(a+k) 


= constant. 


(3.69) 


B(-) 


a 

(a+k) 


(3.70) 


and 


f(-) 


2[n 


(a+kj 

a 


+ 4 1 r| ] = constant. 


(3.71) 


For the case of soft failures, one can derive similar closed-form 
expressions. For example, soft failures with a. ramp reference input r 


J th< T > 


n + 4 | r | 


Ji 


2 2 3 2 

-2(a+k)T if. __aj_ 

5 e 2 4 

2(a+k) (a+k) 3 (a+k) 


2 a 


(a+k)' 


22 2 -(a+k)x 

j-(a+k) T - ■ S . T _ 2a ( (a+k)r+l ) e 

(a+k) J (a+k) 


1/2 


2 ( a+k ) 


5 


(3.72) 
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a 2 <-V 


B(t) 


t) = I i 


c 2 -2(a+k) (i-t f ) 


I T L (a+k) 2 3 

.,2 


2 (a+k) 


(a+k) 


f - . aK . ..2 2K 2 -(a+k)(r-t f ) 

4 (T " t f ) + 77^ (T ' t f ) + « 6 


(a+k)' 


(a+k) 


+ 


2aK 

(a+k) 5 


-(a+k) (-r-t ) 

e r ( (a+k) (x-t f ) 


+ 1 ) 


3K 2 2aK 

(2(a+k) 5 (a+k) 5 


1/2 


(3.73) 


Figure 3.3 shows the results for the hard failure with two different 

estimator speeds. The parameters being used are a=l , t^=.001, c=.01, 

5= . 05 , y=l , K=2 (slow estimator), and K=10 (fast estimator). This shows 
that the threshold selector Is also useful in evaluating the filter dynamics. 
If the dynamics of the filter are fast, there Is a sharper minimum in f and 
generally lower thresholds than when the filter is slow. 


The results for soft failures with the same parameters are shown in 
Figure 3.4. Note that the behavior Is quite different from before. For hard 
failures, there is a single sharp minimum in f which corresponds to a very 
small detection window. For soft failures, f has a hyperbolic type 
behavior. There Is not a unique minimum, and a larger detection window 
compared to the hard failure case needs to be selected. These results are 
quite reasonable as they agree with intuition. 


3.1.3 A Computer-Aided Design Approach for Computing Optimal Thresholds 

It was anticipated that a closed-form solution for the multivariable case 
was not possible hence a computer-aided design (CAD) approach was also 
developed for the threshold selector analysis. The basic operations involved 
In the computation of the thresholds are calculations of frequency response 
characterl sties of L and transient response matrices of L and LT dynamic 
systems. Note that the expression for L is 
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F(tau) 
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Figure 3.3a Minimum Detectable Failure vs. Detection Window (t) 
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Figure 3.3b Minimum Detectable Failure vs. Detection Window (x) 


71 





F(tau) 


SISO SOFT 


“T 

1 












■ 










— 

■ 










■ 










H 

■ 




STEP r 
RAMP r 





! 


n 










5 


■ 









J 

m 

u 








i 


1 

BB 











— — 





- - 


- - - - 








0.0 0.1 0.2 0 . 3 0.4 0.5 0.6 0-7 0.8 0.9 1.0 


Length of detection window (t), sec 


Figure 3.4a Minimum Detectable Failure vs. Oetection Window (t) 
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Figure 3.4b Minimum Detectable Failure vs. Detection Window (r) 


73 




L(s) = (I + F( s) ) 


which is simply the inverse return difference matrix of the filter. To 
compute the transient response of systems with transfer functions l(s) and 
L(s)T(s) , a (minimal) state-space realization of these systems is required. 

In the single-input/single-output (SISO) case, a minimal state-space 

realization for l(s) denoted by ( A ri > 8 rl > c rl > D rl ) can be 

written down by inspection. Similarly, IT is a dynamic system whose transfer 

function involves products of certain inverse return matrices, etc. as shown 

In Figure 3.5. To obtain a minimal realization (A r , B r> C r> 0 r ) for 

this system, the block diagram manipulation facility of CTRL-C [21] called 

INTERC was used. To use this tool, the different blocks in the system block 

diagram are numbered as shown In Figure 3.5. The various Interconnections are 

then specified. The procedure yields a non-minimal realization for the system 

of order n^ = 5 in this case. The procedure MINREAL* is then utilized to 

obtain a minimal realization for the system which is of order n^ = 1 and 

agrees with hand calculations. Suppose h(t) is the impulse response of 

l(s)T(s) system corresponding to the state-space realization (A r , B r> 

C , 0 ) and h,(t) is the impulse response of L(s) system 
r r 1 

corresponding to the state-space realization (A , C r 

0 ). Then the threshold is given by 

rl 


2(n + 5 I r| Q t1 (LT)) 
Q t1 ( L( T-t f ) ) 


(3.74) 


T v U 

2(n + 6|r|(l/* [ f h (t 1 )dt 1 ]-[/ h(t 2 )dt 2 ]dt) 1/2 ) 


f = 


o o 


T~t t 


(3.75) 


(; / f t / h 1 « t 4> dt 4i dt > 1 


n 


*MINREAL uses the staircase algorithm to remove uncontrollable and/or 
unobservable modes. 
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where 


f h(*)d(*) = step response of (A , 8 , C , D ) 
/ r r r r 


(3.76) 


and 


f h (•)<!(•) = step response of (A , B , C , D ) 
I 1 ri rl rl rl 

J o 


(3.77) 


A procedure was then written in CTRl-C to calculate the threshold curve for 
Example 3.1. Figure 3.6 shows the plots of the thresholds generated by CTRL-C 
both for the closed form and the CAD technique using the same values of 
parameters. The integration step for the CAD plot is .001 sec. The slight 
discrepancy in the two curves would disappear if smaller integration steps are 
used. Figure 3.7 shows the plots of threshold selector for a slow (a=l,k=l) 
and a fast (a=l,k= 10 ) estimator. 


Example 3.2 : Multivariable Control of a Turbofan Engine 

To illustrate the idea of threshold selector in the multivariable case, a 
model of a turbofan engine and its multivariable control at sea level static 
conditions with a power lever angle (PLA) of 36° was chosen [12]. 


The states of the system are 

X] » Fan Speed, SNFAN (N^) - rpm 

X 2 = Compressor Speed, SNCOM (N 2 ) - rpm 

X 3 = Burner Exit Slow Response Temperature, 

^t41o " ° R 

X 4 = Fan Turbine Inlet Slow Response Temperature, 
T t4.51o - ° R 

The engine controls are: 

U-] = Main Burner Fuel Flow, WFMB - lb/hr 

U 2 = Nozzle Jet Area, Aj - ft? 
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Figure 3.6 Threshold Selector for Example 3.1 


77 




u 3 = 

U 4 ■ 

U5 = 


Fan Guide Vane Angle, FGV - deg 
Compressor Stator Vane Angle, SVA - deg 
Compressor Bleed Flow, BLC - % 


and the engine outputs are 


Yl - 

Fan Speed, SNFAN (N-|) - rpm 


Y 2 = 

Compressor Speed, SNCOM (N 2 ) - 

rpm 

y 3 = 

Burner Pressure, PT4 - psla 


y 4 = 

Augmentor Pressure, PT 6 - psia 


Y 5 = 

Fan Turbine Inlet Temperature, 

FTIT - 


The normalized system matrices for the example are shown In Table 3.1. The 
open loop poles are at -3.1616, -2.8807, -.7036 and -1.0865. There are two 
transmission zeros at -1.7294 and -.6456 and were computed using the algorithm 
in Ref. 22. It is Interesting that one of the transmission zeros of the 
system is at -.6456 which indicates why It Is not possible to move the slow 
temperature pole at -.7036 very far. The multivariable control law is 
proportional plus integral [6,7] 


G ( s ) = C 
c p 


* s C I 


(3.78) 


where the (normalized) proportional gain matrix (Cp) and the integral gain 
matrix (Cl) are as shown in Table 3.2. The (normalized) filter gain matrix 
(K) is as shown in Table 3.2 and corresponds to filter poles at 


-11.2034, -8.0777, -2.0051, -.6817 


The transfer function matrices of interest are given by 
Nominal Plant: P Q (s) = C(sl-A) 1 B + 0 

Controller; G c (s) = C p + j C T 


(3.79) 

(3.80) 
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Table 3.1 


A = 

- 3 . 9180D+00 
- 1 . 8061D-01 
- 1 . 3190D-01 

- 3 . 8191D-01 


B = 

5 . 1991D-01 
3 . 6266D-01 
2 . 8427D-01 
9 . 3743D-01 


C = 


2 . 2043D+01 
0 . OOOOD+OO 
3 . 7700D+00 
8.0543D+00 
-2 . 9070D+00 


D = 


0. OOOOD+OO 
0. OOOOD+OO 
1 . 0036D+00 - 

9 , 7674D-01 - 

7 . 1316D+00 


System Matrices for the Example 


4 . 1886D+00 
- 2 . 1480D+00 
-2 . 4056D-01 
-1.0501D+00 


1. 1942D+00 
1 . 0836D-01 
3 . 3231D-02 
7 . 3072D-02 


0. OOOOD+OO 
2 . 7339D+01 
1 . 0341D+01 
3 . 1436D-01 
7 . 9884D+00 


0. OOOOD+OO 
0. OOOOD+OO 
8 . 2350D-01 
5 . 7450D+00 
5 . 5560D-01 


-4 . 1148D-02 
1 . 5853D-01 

- 6 . 6630D-01 

- 6 . 7400D-02 


2 . 1974D-01 
7 . 2562D-03 
5 . 7770D-03 
1 . 7417D-02 


0. OOOOD+OO 
0. OOOOD+OO 

- 7 . 6298D-03 
-6 . 6634D-02 

- 5 . 1265D-01 


0. OOOOD+OO 
0. OOOOD+OO 
- 1 . 5200D-01 
- 3 . 8500D-01 
1 . 3247D-01 


1 . 2279D-01 
6 . 6994D-04 
2 . 3770D-04 
-2. OOOOD+OO 


-2 .4990D-02 
- 1 . 2133D-02 
5 . 7672D-03 
2 . 0418D-02 


0. OOOOD+OO 
0 . OOOOD+OO 
-4 . 3237D-03 
- 3 . 7135D-02 
2 . 6855D-03 


0 . OOOOD+OO 
0 . OOOOD+OO 
- 5 . 6233D-02 
9 . 5762D-03 
1 . 5533D-01 


1 . 7226D-02 
7 . 2114D-03 
1 . 6319D-03 
1 . 0634D-01 


0. OOOOD+OO 
0. OOOOD+OO 
5 . 8600D-02 
2 . 2963D-02 
4 . 8290D-02 
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Table 3.2 

Proportional Gain Matrix (Cp), Integral Gain Matrix (Cl), 
Filter Gain Matrix (K) 


K 


2 . 3692D-01 
6.8145D-02 
2 . 3950D-02 
2 . 0241D-02 

1 . 9395D-01 
2 . 8495D-01 
3 . 5780D-02 
6 . 0917D-02 

1 . 0177D-02 
5 . 1618D-03 
2 . 3 L72D-04 
- 5 . 0971D-04 

3 . 7510D-04 
3 . 6090D-05 
2. L026D-05 
-4 . 4286D-06 

-1 . 6788D-02 
-9 . 9733D-03 
5 . 0178D-04 
L . 9605D-03 

Cp 





-4 . 6597D-02 

- 2 . 2281D-02 

- 8 . 5667D-02 
9 . 4676D-02 
3 . 0500D-01 

-2 . 2423D-01 
-4 . 5208D-07 
0. OOOOD+OO 
- 2 . 9942D-01 
2. 3320D+00 

- 5 . 3984D-01 
0 . OOOOD+OO 
0. OOOOD+OO 
0 . OOOOD+OO 
7 . 5286D+00 

-2 . 5951D-02 
3 . 5198D-02 
4. 1154D-02 
5.4087D-02 
0. OOOOD+OO 

0 . OOOOD+OO 
0. OOOOD+OO 
0. OOOOD+OO 
0. OOOOD+OO 
0. OOOOD+OO 

Cl 





- 1 . 0062D+00 
-2 . 3333D-01 
0 . OOOOD+OO 
0 . OOOOD+OO 
0. OOOOD+OO 

0. OOOOD+OO 
0. OOOOD+OO 
0. OOOOD+OO 
0. OOOOD+OO 
0. OOOOD+OO 

0 . OOOOD+OO 
0. OOOOD+OO 
0. OOOOD+OO 
0. OOOOD+OO 
0. OOOOD+OO 

-1 . 7247D-02 
2 . 9998D-02 
0 . OOOOD+OO 
0 . OOOOD+OO 
0 . OOOOD+OO 

0. OOOOD+OO 
0. OOOOD+OO 
0 . OOOOD+OO 
0. OOOOD+OO 
0. OOOOD+OO 
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Filter: 


F(s) 


(3.81) 


From (2.26) 


and 


C(sl-A) 


-1 


K 


G(s) * P q (s) G c (s) 

- C C ( s I —A ) ~ 1 B + 0][C p + \ Cj] 
S(s) = [I + G(s) ] 

T(s) = [I + G(s)]" 1 G(s) 

L(s) * [I + C(sI-A)"'k ] _1 
M(s) = [I + C(SI-A)" 1 K]~ 1 F(s) 


(3.82) 

(3.83) 

(3.84) 

(3.85) 

(3.86) 


max o[L(jw)] = 11.9122 (3.87) 

U> 


Again It is relatively simple to obtain a minimal realization of L(s). The 

state-space matrices are denoted by (A.B„,C.,D„) and are as in 

rfc ri rt r! 

Table 3.3. The poles are at 


-11.2034, -8.077, -2.0051, -.6817 
and the transmission zeros are at 


-3.1616, -.7036, -2.8807, -1.9865 

The block diagram of L(s)T(s) is as shown in Figure 3.8. INTERC was 
used to obtain a non-minimal realization of the system. Note that for the 
purposes of defining the output to CTRL-C, it is necessary to introduce a 
fictitious block (6) in the diagram. The order of the non-minimal realization 
is n $ = 27. MINREAl was used to yield a minimal realization of order 
n $ = 10*. The final results for (A^, 8 r< C^, 0^) are displayed in Table 3.4. 
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Table 3.3 


Associated State-Space Matrices (A ri , B r ^, C ra , D,^) 


ARL 

-2 . 0003D+00 
- 7 . 2277D-03 
- 1 . 2263D-01 
-1 . 5242D-03 


BRL 

- 2 . 0366D-02 
2 . 3844D-02 
2. 3692D-01 

6 . 8145D-02 


5.9409D-02 
-6 . 6573D-01 
- 5 . 0294D-02 
1.5345D-01 


-6. 1104D-02 
3 . 5460D-02 
1 .9395D-01 

2 . 8495D-01 


8 . 2386D-01 
-6 . 5511D-01 
- 9 . 2306D+00 
- 1 . 7314D+00 


5 . 0848D-04 
2 . 3439D-04 
1 . 0177D-02 

5 . 1618D-03 


2 . 7009D+00 
-1. 2030D+00 
-1 . 3534D+00 
-1 . 0071D+01 


4 . 3184D-06 
2 . 1049D-05 
3 . 7510D-04 

3 . 6090D-05 


- 1 . 9631D-03 
4 . 9150D-04 
- 1 . 6788D-02 

-9 . 9733D-03 


CRL 

0 . 0000D+00 
2 . 2204D- 16 
4 . 3636D-03 
3.7484D-02 
- 8 . 3267D- 17 


1 . 3878D- 16 
3 . 3307D-16 
-7.6070D-03 
-6 . 6438D-02 
-5 . 1265D-01 


2 . 2043D+01 
2 . 2204D-16 
3 . 7700D+00 
8 . 0543D+00 
- 2 . 9070D+00 


6 . 6613D- 16 
2 . 7339D+01 
1 . 0341D+01 
3 . 1436D-01 
- 7 . 9884D+00 


DRL 


- 1 . 

0 . 

0 . 

0 . 

0 

0 . 

- 1 . 

0 . 

0 . 

0 

0 . 

0 . 

- 1 . 

0 . 

0 

0 , 

0 . 

0 . 

- 1 . 

0 

0 . 

0 . 

0 . 

0 . 

-1 
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The poles of A p are 

-11.2034 + O.OOOOi 
-8.0777 + O.OOOOi 
-4.7046 - O.OOOOi 
-2.6196 + 1 . 8997 i 
2.7863 + O.OOOOi 
-2.0051 + O.OOOOi 
-0.7347 + O.OOOOi 
-0.6817 +■ O.OOOOi 
-0.1034 +■ O.OOOOi 

and the transmission zeros are 


-1.9865, -.7036, -3.1616, -2.8807 


which are as expected. 


The (0 ) LT is then of the form 

t1 

T t t 

(Q t1 )LT(t) = ( [ j h(t 1 )dt 1 ] T [ y h(t 2 )dt 2 ]dt j 1 /2 (3.88) 

0 0 0 

where h(*) is the impulse response matrix of the (A r> B r , C r> 

D ) system and its (i.j)-th element Is the response of the i-th 
component of the output to an Impulse applied at the j-th component of the 
input while all other components of the Input are zero and the Initial state 
is zero. Note that an expression such as: 



V dt ! 


step response matrix 


(3.89) 


* This is a numerically delicate problem. The tolerance for MINREAL should be 


chosen to be 100 


A 8 
C D 


c where c Is the machine precision. 

2 
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Table 3.4 


Final Results for (A r , B r> C r , 0 f ) 


AR 


Starting at row 1 columns 1 thru 


1 .8833D+00 
3 . 0702D-02 
3 . 0137D-01 
1 .6443D-01 
1 . 3275D-01 
3 . 1541D-02 
0 . OOOOD+OO 
0 . OOOOD+OO 
0. 0000D+00 
0. OOOOD+OO 


2 . 2461D-02 
- 1 . 9773D+00 
1 . 1797D-02 
- 2 . 1904D-01 
8 . 5894D-03 
- 1 . 1947D-02 
1 . 2095D-01 
0. OOOOD+OO 
0. OOOOD+OO 
0. OOOOD+OO 


4. 1244D-01 
5. 3534D-02 
- 7 . 7122D-01 
7 . 2472D-02 
-5. 1391D-02 
-3 . 1920D-02 
- 1 . 6299D-02 
1 . 4931D-01 
0. OOOOD+OO 
0. OOOOD+OO 


6 

-9 . 7421D-02 
- 1 . 2298D-01 
2. 5823D-02 
- 6 . 8127D-01 
3 . 4825D-04 
-2 . 5310D-02 
-4 . 6428D-02 
1 . 8504D-02 
1 . 5650D-01 
0 . OOOOD+OO 


8 . 6466D-02 
-1.4802D-01 
8 . 9131D-01 
1 . 8765D+00 
4 . 0643D-01 
2 . 0529D+00 
3 . 1668D-01 
1 . 4743D+00 
-7 . 9998D-01 
1 . 1380D+01 


Starting at row 1 columns 7 thru 


4 . 7572D-01 
6.1265D-01 
-2 . 1959D+00 
-2 . 7577D+00 
-1. 1208D+01 
-4 . 9271D+00 
-9.4908D+00 
-4.2548D+00 
6. 6487D-01 
- 1 . 2398D+01 


1 . 7422D+00 
-1 . 5215D+00 
- 1 , 3057D+00 
6 . 3804D-Q1 
1 . 9S45D+00 
1 . 3970D+00 
4 . 9505D-02 

- 8 . 4370D+00 

- 6 . 2387D-01 
3 . 6584D+00 


2 . 1100D+00 

- 9 . 1115D-01 

- 2 . 2261D+00 

- 2 . 9244D+00 

- 6 . 1724D+00 
- 1 . 6149D+00 

- 6 . 5492D+00 
- 1 . 4728D+01 

- 6 . 3950D+00 
- 1 . 6509D-01 


10 

4 . 8148D-01 
1 . 1621D+00 

- 5 . 5108D-01 
-7 . 6130D-01 
- 1 . 0306D+00 
- 1 . 3658D+00 

- 3 . 6488D-01 

- 6 . 2005D-01 
1 . 2374D-01 

-6 . 3510D+00 


BR 


1 . 2913D-02 
3 . 5216D-02 
-6 . 2897D-02 
* 1 . 1748D-01 
-2 . 6848D-01 
- 1 . 8192D-01 
-5. 3986D-02 
- 1 . 3163D-01 
5.2728D-02 
- 6 . 2446D- 01 


-3. 5065D-03 

- 3 . 5419D-03 
5 . 7229D-03 

-2 . 7295D-02 
1 . 7852D-02 
-1.4634D-03 
-4 . 5738D-02 

- 5 . 5903D-04 

- 5 . 1607D-02 

- 7 . 4094D-03 


6 . 7581D-02 
9 . 7594D-02 

- 9 . 2170D-03 

- 2 . 5362D-02 
5 . 2850D-02 

-1.2483D-03 
-1 . 1493D-01 
6 . 5713D-03 
- 1 . 2889D- 01 

- 2 . 7436D-02 


- 9 . 3985D-03 
- 1 . 3254D-02 

1 . 8487D-03 

- 5 . 6664D-03 
-2 . 0758D-02 

1 . 4683D-02 
2 . 1870D-02 
-4 . 2030D-03 
8 . 8618D-03 
1 . 8488D-04 


0. OOOOD+OO 
0. OOOOD+OO 
0. OOOOD+OO 
0. OOOOD+OO 
0. OOOOD+OO 
0 . OOOOD+OO 
0 . OOOOD+OO 
0 . OOOOD+OO 
0. OOOOD+OO 
0 . OOOOD+OO 


6 . 8098D-04 
4 . 3156D-02 
1 . 5985D-01 
3 . 0526D-01 
6 . 3688D-01 
4.4477D-02 
7 . 7116D-01 
1 . 9028D-01 
1 . 4590D-01 
2 . 1815D+00 
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Table 3.4 (Continued) 


CR 


row 1 co 

- 3 . 3858D- 16 

- 3 . 4001D- 16 
-9 . 3675D- 17 

- 3 . 1919D- 16 
2 . 7756D- 16 

Starting at row 1 co 
-3 . 7131D+01 - 3 . 0531D- 16 

9 . 5083D+00 - 2 . 3928D+01 

-2 . 7291D+00 -1.0432D+01 

- 1 . 5132D+01 -5.4949D-01 

- 3 . 3307D- 16 2 . 3384D- 15 


lumns 1 thru 6 

-1.1796D-16 -3 . 6082D-16 

-1.9663D-15 -1 . 3704D-16 

- 2 . 6877D* 16 -1.7781D-17 

-2 . 2291D-16 -I . 7347D-16 

2 . 6472D-15 -4.7184D-16 

lumns 7 thru 

- 8 . 8840D+00 

- 3 . 9740D+01 
-2 . 3954D+00 

4 . 9940D+00 
7 . 6032D+01 


Starting at 
1 . 9429D- 16 
- 5 . 0220D- 16 
-6.8955D-17 
1.4572D-16 
1 . 3878D- 16 


0 . OOOOD+OO 
0 . 0000D+00 
3 . 7478D-02 
- 9 . 9875D-02 
1. 8695D-01 


0. 0000D+00 
0. OOOOD+OO 
1 . 6983D-01 
1 . 2594D-01 
9 . 2186D-01 


0. OOOOD+OO 
0. OOOOD+OO 
4 . 8632D-01 
2 . 8937D-01 
1 . 7470D+00 


10 

0 . 0000D+00 
1.5959D-16 
1 . 0780D+00 
- 2 . 7756D-17 
-6 . 9389D-17 


0 . OOOOD+OO 
0 . OOOOD+OO 
2 . 6586D-02 
1 . 8044D-01 
3 . 1605D-02 


4 . 4409D- 16 
-9 . 0206D- 16 
-5 . 8981D- 17 
1 . 6653D- 16 
6 . 1062D- 16 


0. OOOOD+OO 
0 . OOOOD+OO 
0 . OOOOD+OO 
0. OOOOD+OO 
0 . OOOOD+OO 


2 . 2204D-16 
6 . 9389D-18 
1 . 0446D+00 
6 . 8199D+00 
1 . 6653D- 16 
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Is a step response matrix, where the (l.j)-th element of it is the response 
at time t of the i-th component of the output when the j-th component of 
the input is a step function while all other components of the input are zero 
and the initial state Is the zero state. In this case the step response 
matrix Is a 5x5 matrix. Similar statements apply to (Q^JL operator. 

We finally compute the multivariable threshold as 

2J th ( T ) 

f (T)= B(TT (3 - 90) 


where 


J th (r) = (max o[L(j«)])n + 6 | r | oC ( 0 xl LT) { x) ] (3.91) 

td 

B(t) = o[(O tl l)(T - t f )] (3.92) 

linn < n (3.93) 

®(A(j«)) < 6 (3.94) 


The threshold selector results for the hard failure case are as shown In 
Figure 3.9. Note the similarity with the scalar case. The results for the 
soft failure case are as shown in Figure 3.10. 

3.1.4 Summary 

This chapter has discussed the effect of model error on detection. An 
innovative framework was developed to determine the effect of model error on 
sensor failure OIA algorithms analytically. A new concept called threshold 
selector was introduced. The threshold selector analysis allows determination 
of optimal threshold and size of smallest detectable failure as a function of 
model error bound, noise, variance, speed of the filter, class of reference 
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Length of detection window t (sec) 


Figure 3.10a Threshold - Soft Failure Case, Example 3.2 
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Input signal, and class of failure signal. The threshold selector constitutes 
a powerful design aid and allows one to arrive at a proper balance between 
robustness and OIA sensitivity. The analysis In this chapter can be extended 
to Isolation and accommodation problems, as well as other DIA designs than 
those discussed here. 



IV. ROBUST FILTER DESIGN FOR DIA 


This chapter discusses the design of robust filters for use as part of a 
sensor failure DIA algorithm. The robustness of the filter Is due to 
Introduction of an internal model as well as frequency-shaping of the LQG cost 
functional. The necessity of the internal model was discussed in Chapter II 
(see Section 2.2.3). Recall that the presence of an internal model results in 
asymptotically unbiased filter estimates. The Internal model also provides 
robustness with respect to parameter perturbations. However, since other 
sources of model uncertainty are present (see Chapters I and II), the filter 
may be made more robust by frequency shaping of the LQG cost functional. This 
Is a result of taking the bound on model error Into account. The concepts are 
applied to a multivariable turbofan engine example. 

4.1 PROBLEM FORMULATION 

A filter may be made robust by introduction of an Internal model and by 
frequency shaping . The internal model provides robustness with respect to 
parameter perturbations and results in asymptotically unbiased estimates. The 
filter may be made more robust by adding dynamics to the filter to compensate 
for other types of modeling uncertainty. This can be done in a formal way by 
replacing the constant weighting matrices in a standard LQG cost functional 
with weighting matrices which are functions of frequency. This is referred to 
as frequency shaping. The weighting matrices are chosen to reflect model 
uncertainty. For example, if there is unmodeled high-frequency dynamics, the 
weighting matrices may be chosen to be constant over the frequency range where 
the system model Is known accurately and increase as a function of frequency 
In the frequency range where the model is less accurate. This section 
discusses the design of robust filters for sensor failure DIA which employ 
both an Internal model and frequency shaping. 
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4.1.1 Filter with Internal Model 


The Internal model principle and its application to robust filter design 
was discussed in Section 2.2.3. It was shown that the internal model is 
necessary to provide asymptotically unbiased estimates and robustness with 
respect to parameter perturbations. For a system with unknown constant 
measurement biases, an Internal model (Integrators) may be Introduced by state 
augmentation: 
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(4.1) 


(4.2) 


where KP and KI are the proportional and integral estimator gains 


4.1.2 Filter Design with Frequency Shaping 


Consider the system described by: 


state: x * Ax + Bu +• B-|W (4.3) 

measurement: z = Cx + Du + v (4.4) 

where w and v are independent, zero-mean, white Gaussian process and 
measurement noise processes. The standard Kalman filter for this system Is of 
the form 

x = Ax i- Bu + K(z - Du - Cx) (4,5) 

where K Is the Kalman filter gain. We desire to find a filter which 
minimizes the performance index. 
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(4.6) 


GO 

$ m \ J {(z - Cx) T R~\z - Cx) + w T 0 _1 w} dt 


*2 J ( w ^ 0 ] w + v T 


R 1 v} dt 


In this cost function, Q and R are constant matrices Independent of 
frequency. Using Parseval's theorem, the performance index In (4.4) can be 
transferred to the frequency domain 


*/ “ ^ f (W*(ju) Q Vjw) + V*(J«) R _1 V(jw)} du 


(4.7) 


Note that the two terms in the above integrand have the constant weighting at 
all frequencies. However, the model may be well known within a certain 
frequency range and not known accurately outside that frequency range. It 
would then seem desirable to have weighting matrices which are functions of 
frequency and be able to reduce the filter gain outside model bandwidth to 
reduce sensitivity and increase performance of the filter. It is possible to 
consider making Q and R functions of frequency 

/ CO 

{W*( ju) Q -1 (u 2 ) W (ju) + V*(ju) R -1 (w 2 ) V( jw) } du (4.8) 

- CD 

A sufficient condition for the existence of a filter minimizing (4.8) is 

2 

that Q(ju) and R(ju) be positive semi-definite matrices in J . 


The problem as posed in (4.8) can be converted to a standard LOG problem 
as shown below. One can treat w and v as colored noise sources generated 
by shaping filters of the form 


W(ju) = Q 1/2 (ju) w'(ju) 


V(ju) = R 1 / 2 (ju) v'(ju) 


(4.9) 

(4.10) 


where W' and V' are white noise processes. Q 1/2 and R 1/2 are square 
roots of 0 and R such that 
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(4.11) 


1/2 1/2 * 

Q(J«) = 0 (j«) CO (jo,)] 

R( jo») = R 1/2 (j<o) [R 1/2 (jo,)]* (4.12) 

i /? i/2 

where Q (jo>) and R (jo>) are rational functions of jo,. 


In sensor failure 0 I A . the primary cause for the use of frequency-shaping 
Is unmodeled dynamics. In our application to the engine problem, Q is taken 
to be constant ( 1 . e . , Independent of frequency) and R Is chosen as 


, 2 

r 2 (o> + r^) 

r l (w 2 + r 2 ) 


> r. 


(4.13) 


(see Example 4.1). Note that strictly proper transfer functions in 
-1 /2 

R (j«) would cause R(jo>) to approach zero at high frequency. This 

implies perfect measurements. Therefore, in practice, one should choose 

- 1 / 2 . . . . 

R (ju) to be proper. 


Next, the modified measurement zi is considered 
Z] = z-0u = Cx + v = y + v 
If we define a shaped-measurement vector 

Z-j ( j<«>) = R(jwf 1/2 Z^jw) 


(4.14) 


(4.15) 


and the noise-free shaped output 


Y ' ( jw) = R(jo,f 1/2 Y ( jw) 


(4.16) 


-1 /2 

then a realization of the system with transfer function matrix R(jw) is 


x = A x + B„ Cx 
v v v v 


(4.17) 
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y = c v x v + ° v Cx 

(4.18) 

with the shaped-measurement equation 


Z 1 = C v x v + D v Cx + v 

(4.19) 

where v' Is a white noise process as in (4.10). 


The combination of equations (4.3) and (4.4), and (4.17), (4.18), and 
(4.19) defines a dynamic system driven by independent white noise sources and 
is a standard Kalman filtering problem: 

• <1 
x = Ax + Bu + K e (z ] - z 1 ) 

(4.20) 

• 

x v = A v x v + 8 v Cx + K v (z l ' V 

(4.21) 

z 1 is obtained from z using the realization of 

R" 1/2 : 

*z • \ x z* B v Z 1 

(4.22) 

i 

Z, = C X +0 Z, 
1 V Z V 1 

(4.22) 

Since equations (4.17) and (4.18), and (4.22) and (4.23) involve the same 
realization, the redundancy in the estimation equations can be eliminated by 


defining a new set of states [23] 

x ‘ " x z ' x v (4.24) 

The overall frequency-shaped estimator equations are then given by 
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(4.25) 


x = [C : 0] 


(4.26) 


The control law is then based on x 


u - K R | x ! 


(4.27) 


where K R Is the regulator gain. A block diagram of the system showing the 
frequency-shaped optimal estimator is shown in Figure 4.1. The filter in Eqs 
(4.22) and (4.23) acts as a prefilter on the measurements. 

-I /2 

Theorem 4.1 : The zeros of R are the transmission zeros of the 

frequency-shaped estimator whenever the sensors are frequency-shaped 
Individually. 


Proof : The transmission zeros of the estimator (with input z and 

output x) are given by those frequencies such that the matrix pencil 
S(x) loses rank 
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< X 



TOO 


Figure 4.1 Block Diagram of Estimator (Eqs. 4.25 and 4.26) 
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S(X) = 

- (K v ° v - B y ) C 

S' - \ 4 \ 
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I 

0 

0 


(4.28) 


(4.29) 


By performing elementary row and column operations, S(x) is equivalent to 


S(X) 


Si - A j 
1 

0 j (si 

i 

I 

1 i 


K D 
.1 


A +8 D 1 C ) 
v v v v' 


e v - "v 


(4.30) 


which loses rank at the location of the eigenvalues of (A y + B y D y C y ) 

- 1/2 

which are the zeros of R since it is diagonal, as are A y . B v , C y , 

and D . Q.E.D 
v 

- 1/2 

In the case of individually frequency-shaped sensors, R would be 
diagonal and, if a frequency shaping of first order is introduced in all 
sensors then 


-1/2 P^s + 

R ii ° Z^TTTjT 


(4.31) 


and hence 


A v ■- diag (P.) 


B - I 
v 


C y = diag 


P 1 (Z 1 ~ P j> 


(4.32) 

(4.33) 


(4.34) 
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(4.35) 



Notice that R ^ ^ ' s have unity D.C. gain. These will Introduce lag 
compensators into the measurement loops provided that 

IZ. | > IP^ (4.36) 

which implies more noise at higher frequencies. Frequency shaping then 
improves the robustness properties of the system at higher frequencies, while 
maintaining performance at low frequencies. Notice that this is similar to 
introduction of first-order lags of classical control, but is inherently 
multivariable. 

4.2 EXAMPLE 4.1 

This example illustrates the ideas of internal model and frequency 
shaping presented in Sections 2.2.3 and 4.1, respectively. The internal model 
provides for asymptotically unbiased estimates in the presence of biases and 
parameter errors, l.e., the estimate of the engine outputs track the output 

measurements (i.e., z - z -» 0). Frequency shaping provides for robustness 
with respect to unmodeled dynamic uncertainties. The combination of the 
internal model and frequency shaping results in the desirable robustness 
properties of the filter. 

A turbofan engine model and its multivariable control law (the same as in 
Example 3.2) at sea level static condition and PLA=36° was chosen for design 
purposes. Figure 4.2a shows a steady-state run corresponding to Revision 2 of 
the previous program [9]. This figure shows that the estimate of N1 , N2, PT4 , 

PT6, and FTIT do not track the measurements (i.e., z - z * 0 asymptotically). 
If an internal model (integrators) for biases is introduced, then the 
augmented system has the form of (4.1) and (4.2). The estimator gains were 
computed using CTRL-C [21] and are shown in Table 4.1. Figure 4.2b shows a 
steady-state run with the internal model present in the engine/control model. 
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Figure 4.2a Engine Steady-State Response at Sea Level Static Condition and 
PLA * 36° with No Internal Model Present (Revision 2, Ref. 9) 
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Figure 4.2b Engine Steady-State Response at Sea Level Static Condition and 
PLA = 36° with Internal Model 








Table 4.1 


Kalman Filter Gain Matrix 


KK 


8. 3760 

3. 1891 

0 

4. 5611 

12. 4951 

0 

1. 3178 

2. 0933 

0 

0. 0701 

-0 0033 

0 
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l . 0903 

0 

0. 0009 

2. 9770 
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0. 0121 

0. 0476 

1 

0. 0019 

0 0036 

0 
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-0. 1824 

-0 


2311 

0. 0303 

-0. 6314 

4499 

0. 0173 

-1. 2775 

0836 

0. 0048 

-0. 2963 

0006 

0. 0001 

-0. 0004 

0253 -O'. 0018 

-0. 0429 

OP 5 1 
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0. 0000 
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-0 0003 
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-0. 0001 

1. 0077 
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(This corresponds to Revision 3 of the previous program [9].) Figure 4.2b 
shows that the estimates now asymptotically track the measurements (i.e., 

z - z -* 0) for all outputs. Note that there is a dip in Figures 4.2a and 
4.2b, which is due to initialization transients of the engine simulation. 
Figures 4.3a and 4.3b are transient runs at sea level static condition and a 
PLA step from PLA=36° to 52°. Figure 4.3a shows transient responses with no 
Internal model and Figure 4.3b illustrates the responses with the internal 
model present. It can be seen from Figure 4.3b that the estimates follow the 

measurements, I.e., z -» z ) 0, in steady state. 

There are now two Internal models in the closed-loop system (i.e.. 
Integrators both In the filter and controller. Figure 4,4 is a block diagram 
of the system showing the two Internal models. The presence of the internal 

model in the filter ensures z - z -* 0, and the presence of the internal 
model in the controller ensures z - r -» 0, which implies that r - z -* 0 in 
steady state. Note that the controller has a partial internal model, I.e. it 
only has Integrators on N1 and PT6 outputs (see form of Cl in Table 4.2). 
Therefore, even though all the estimates are unbiased, we can only guarantee 
zero steady-state tracking error in N and PT as seen in Figure 4.2b 

and 4.3b. Note that in Figure 4.2b and 4.3b for the PT6 output, z - z error 
has become zero whereas z - r error has not. This is because the Integrator 
(in the control law) associated with PT6 output has a small gain which 
explains why this error is slow in decaying to zero (for details, see Ref. 9). 


We can now proceed to add frequency shaping in the filter. Based on the 
results of a bound on model error (Appendix A), we chose to frequency-shape 
all sensors using 


-1/2 10ir(s + 20? ) 

K 1 i {5) 20ir( s + 10*) 


(4.37) 


I.e., a first-order lag with breakpoint at 5 Hz. The frequency-shaped system 
matrices ( A v , B^ C^, D^} are shown in Table 4.3, and Table 4.4 shows the 
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Engine Transient Response at Sea Level Static Condition and PLA 
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Model Present (Revision 2, Ref. 9) 
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Figure 4.3b Engine Transient Response at Sea Level Static Condition and PLA 
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Figure 4.4 Overall Robust Servomechanism with Internal Models 



Table 4.2 


Proportional and Integral Gain Matrices 
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Table 4.3 

Frequency-Shaped System Matrices 
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Table 4.4 

Frequency-Shaped Filter Gain 
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frequency-shaped Kalman filter gain matrix. The closed-loop filter poles are 
shown In Table 4.5. Figure 4.5 shows the transient engine response corres- 
ponding to the closed-loop system with both internal models and frequency 
shaping. Figure 4.5 shows the effects of the frequency shaping in this case, 
i.e., the engine transient response has been slowed down slightly. 

4.3 SUMMARY 

In this chapter, we have discussed the design of robust filters for OIA. 
The robustness properties of the filters are twofold. First, the filter is 
made robust with respect to parameter perturbations, using an Internal model. 
This Is an extension to internal model principle of multivariable robust 
servomechanism theory. The robustness property achieved Is due to creation of 
certain structurally robust blocking zeros. Second, the filter was made 
robust with respect to other sources of uncertainty via frequency shaping. 

This robustness property is achieved also through creation of certain 
transmission zeros. The results were applied to a multivariable turbofan 
engine example. 


112 



Table 4.5 


Closed-Loop Filter Poles 


-23. 4786 

+ 

4. 921 5 i 
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0. 0000 i 

-1. 8294 

- 
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Figure 4.5 Transient Response at Sea Level Static Condition and PLA Step 
Input from 36° to 52° at 0.2 Seconds, and with Internal Model 
and Frequency Shaping 
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V. EVALUATION RESULTS 


5.1 INTRODUCTION 

This chapter evaluates the results of this program through the validation 
of the threshold selector results of Chapter III. Before the actual 
validation. It Is useful to summarize the results of the report and put the 
evaluation results into perspective. 


This report has presented the results of recent research in the 
development of robust fault detection, isolation, and accommodation ( D I A ) 
algorithms for sensor failures. Specifically, tools and procedures have been 
developed that allow a designer to use Information about model uncertainty 
when designing a sensor DIA logic. This is a major step In developing an 
ability to design and Implement a practical fault-tolerant control system. 

A DIA system, as treated In this report, consists of three main 
components : 

(1) a filter that compares measurements to predictions (based on a 
model) to produce an innovations sequence; 

(2) a norm computation that reduces the innovations to a single 
measure useful for comparing against a threshold; and 

(3) a threshold. 

The goal of the design process described in Figure 1.9 is to select a 
combination of these three components to produce a system that has adequate 
performance (smallest magnitude of failure detectable, speed of detection, and 
minimum false alarm rate) without being excessively complex. 

The emphasis of this reported effort has been to provide tools and 
procedures that allow the design process of Figure 1.9 to be carried out. 
Specifically produced have been techniques that 
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(1) permit the performance of filter innovations measure combination to 
be computed analytically using an estimate of the model uncertainty 
in the system; and 

(2) permit a filter's performance to be improved by Incorporating a 
knowledge of the modeling error bounds in its design. 

The first was addressed in Chapter III under the heading of the Threshold 
Selector; the second was addressed In Chapter IV. 

Selecting the best f llter/innovations measure combination is a multistep 
and Iterative process (Figure 1.9). The filters used to generate the 
innovations sequences can vary in complexity from constant gain types to ones 
that Include frequency shaping and Internal models to a fully adaptive design 

As a general rule, increased complexity Is required to improve the filte 
performance. The study reported in Ref. 9 dealt with constant gain filters. 
This report described the use of frequency shaping and internal models to 
improve filter performance in the presence of modeling errors. Adaptive 
filter designs are left for future studies. 

Many measures of the size of innovations sequence are also possible for 
consideration. Examples include weighted sum squared residual (WSSR), 
likelihood ratio, and generalized likelihood ratio. The performance of each 
in combination with the different filters could be different and should be 
investigated before a final system Is designed. Note, however, that the only 
measure of the innovations sequence dealt with In this effort in the WSSR 
norm. This is the same as that used in Ref. 9. 

While the filter and size of innovations measure selections are clearly 
critical to designing a successful OIA system, it Is the ability to evaluate 
analytically the performance of the combination that makes the design process 
of Figure 1.9 feasible. This capability is what is provided by the Threshold 
Selector described in Chapter III. It produces an estimate of the sizes of 
the smallest failures that can be detected and a measure (i.e., threshold) 
against which the norm of the Innovations sequence can be compared to 
determine the presence of a failure. 
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As an example, the threshold selector can be used to predict the 
performance of the robust filter designed in Chapter IV. Figure 5.1 shows the 
threshold selector results for a soft failure using the robust filter with 
Internal model and frequency shaping designed in Example 4.1. This figure 
shows that the robust filter is capable of detecting failures of smaller size 
compared to the constant gain filter developed in Ref. 9 (compare Figures 5.1 

and 3.10) . 

5.2 VALIDATION OF THE THRESHOLD SELECTOR 

Because of its critical role in the design process, the ability of the 
Threshold Selector to predict realistic thresholds has been validated 
experimentally. The goal was to demonstrate the fact that induced failures 
are detected and that false alarms are avoided. 

The filter/norm chosen for this demonstration Is the same as that 
developed in Ref. 9 for a multivariable turbofan engine. Specifically, the 
filters are constant gain filters with no frequency shaping or internal 
models. The norm is a WSSR. The flight condition is sea level static at 36° 

PLA. 


The reason for this filter/norm choice is that it allows direct 
evaluation of the Threshold Selector only. Reasonable thresholds were 
determined empirically for this combination and a full evaluation was 
performed in a previous program [9]. Consequently, there is a data base 
against which to compare the results obtained with a new threshold Implemented. 

Compared In Table 5.1 are the thresholds determined empirically in the 
previous program [9] and the thresholds computed with the Threshold Selector. 
Note that while the results are of similar magnitudes, the Threshold Selector 
computed values are smaller. This indicates that failures of smaller 
magnitude can be detected (and faster) with the filter/norm combination than 
were previously expected. Required to be validated experimentally Is that 
false alarms are not induced as a result of decreasing the thresholds. 
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Table 5.1 


Comparison of Sizes of Thresholds 


TYPE OF FAILURE 

SIZE OF THRESHOLD USED IN 
REVISION 2 OF REF. 9 

SIZE OF THRESHOLD COMPUTED 
BY THRESHOLD SELECTOR 

HARO 

2.0 

.62 

SOFT 

1 .43 

.62 


Our evaluation consisted of the comparison of failure detection times for 
various failures using the thresholds in Ref. 9 and the thresholds computed 
using the threshold selector. A turbofan engine dynamic model (the same one 
as In Examples 3.2 and 4.1) at sea level static condition and PLA = 36° was 
used for this evaluation. Figures 5.2 through 5.4 compare the response of the 
same DIA algorithm (i.e., the same filter/norm combination) for the two 
threshold levels presented in Table 5.1. Figures 5.2 and 5.3 present selected 
responses to a hard failure in N1 for the 0 I A scheme using the empirical 
threshold and the Threshold Selector computed threshold, respectively. 

Figures 5.4 and 5.5 present the results for a hard failure in N2. Figures 5.6 
and 5.7 present the results of a drift failure In N1 . Figures 5.8 and 5.9 
present the results of a drift failure in PT4. Note that key events 
characteristic of all the plots are Indicated in Figure 5.2. 

The results are as expected and are summarized In Table 5.2. Large step 
failures create a WSSR norm of the innovations sequence larger than the 
threshold in both cases and trigger a failure Indication within one window 
width of the WSSR norm. The size of the norm has no effect on the 
performance. For drift failures, however, the smaller threshold does permit 
more rapid detection of the failure. 
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Figure 5.2 (Continued) 
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Figure 5.3 (Continued) 
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Figure 5.4 (Continued) 
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Figure 5.5 (Continued) 
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igure 5.6 Failure Transients for N1 Sensor Drift Rate of 300 RPM/second at 
Sea Level Static Condition and PLA = 36° using Method [9] (see 
Table 5.2) 
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Figure 5.7 Failure Transients for N1 Sensor Drift Rate of 300 RPM/sec at Sea 
Level Static Condition and PLA = 36 using Modified Thresholds 
(Table 5.1) Computed from the Threshold Selector (see Table 5.2) 
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Figure 5.8 failure Transients for PT4 Sensor Drift Rate of 10 PSI/sec at Sea 
Level Static Condition and PLA = 36 using Method [9] (see Table 
5 . 2 ) 
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Figure 5.9 Failure Transients for PT4 Sensor Drift Rate of 10 PSI/sec at Sea 
Level Static Condition and PLA = 36 using Modified Thresholds 
(Table 5.1) Computed from the Threshold Selector (see Table 5.2) 









Table 5.2 


Comparison of Detection Times 


TYPE OF FAILURE 

DETECTION TIME (SEC) USING 
METHOO (REVISION 2) [9] 

DETECTION TIME (SEC) 
USING THRESHOLD SELECTOR 

N1 HARO FAILURE 
+1000 RPM STEP 

0.002 

0.002 

N2 HARD FAILURE 
+1000 RPM STEP 

0.002 

0.002 

N1 DRIFT RATE OF 
+300 RPM/SEC 

4.466 

1 .734 

PT4 DRIFT RATE 
OF +10 PSI/SEC 

4.584 

2.042 


Of primary importance in this comparison is that this increase in 
performance was gained without a penalty of false alarms. This verifies that 
the Threshold Selector does produce a useful estimate of the threshold and can 
replace the previously required empirical approach to threshold calculations 
to accommodate modeling errors in addition to noise considerations. 

5.3 SUMMARY 

This chapter has presented evaluation results for this research effort. 
The results of this study were summarized at the beginning of this chapter. 

The evaluation validated the threshold selector results. It was verified that 
the generally lower thresholds, predicted by the threshold selector analysis, 
results in improvements for soft failure detection. This is done without 
triggering false alarms. The detection time was lowered by at least a factor 
of two for soft failures compared to previous techniques. 
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VI. CONCLUSIONS 


This report contains the results of recent research in the area of robust 
fault detection, isolation, and accommodation for sensor failures. The 
results have been illustrated by application to sensor DIA for an aircraft 
turbine engine. At the center of attention Is model uncertainty. Model 
uncertainty has been singled out as the main source of problems with previous 
DIA algorithms as it affects performance. Various sources of model 
uncertainty were discussed. The effects of model uncertainty are represented 
by a bound as a function of frequency. This Is consistent with recent 
approaches In robust control theory. The effects of model uncertainty on 
stability, asymptotic tracking, and tracking performance was studied and new 
performance robustness measures were derived. The effects of model 
uncertainty on failure detection were also studied. The main machinery used 
was robust multivariable control theory. Fundamental results were derived for 
selection of optimal thresholds in innovations-based DIA algorithms. The 
estimator logic used in the DIA technique was made robust by providing 
robustness both at low and high frequencies. Evaluation results show 
improvements compared to previous techniques. 

The general contribution of this research has been the extension of 
recent advances in robust control system design to sensor DIA and estimator 
design. The specific contributions are: 

- analysis tools with which to quantify the trade-off between 
performance robustness and DIA sensitivity; 

- design methods which allow higher levels of performance 
robustness to be achieved for given levels of DIA sensitivity; 

- demonstration of the applicability of these tools using an 
aircraft turbine jet engine multivariable control example. 

A powerful synthesis tool has been developed for DIA algorithms. This would 
allow for optimal achievable levels of performance. In particular: 
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- a "threshold selector" has been created which quantifies the 
effects of uncertainty on DIA performance; 

- measures have been derived to quantify the uncertainty and the 
performance robustness. 

The advances to robust estimator design to achieve higher levels of 
performance robustness include 

- the development of estimators using the "internal model 
principle" to achieve asymptotic convergence despite model error; 

- the incorporation of frequency weighting in an LOG cost 
functional to modify an estimator design for DIA. 

The results were demonstrated on a dynamic simulation of a multivariable 
turbofan jet engine example and showed improvements over previous techniques. 

The results of this research can be pursued in other directions. It is 
possible that even higher levels of performance are achievable in some cases 
by an adaptive technique (see Figure 1.9). It would also be Interesting to 
pursue the idea of this research in a decentralized control strategy. This Is 
Important as integration of flight and propulsion control systems is being 
considered [25]. 
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APPENDIX A 

GENERATION OF BOUND ON MODEL UNCERTAINTY 


A. 1 SYSTEM UNCERTAINTIES 

All nominal design models of a system (plant) contain some degree of 
modeling errors. High fidelity models represent a plant more accurately than 
others. These errors are called "model uncertainties." There are, In 
general, two types of uncertainties: structured and unstructured. For 
example, the former refers to model parameter uncertainties, whereas unmodeled 
dynamics lies in the latter category. Reduced order modeling, linearization 
about operating points, neglecting nonlinearities, all result in contributions 
to either structured or unstructured uncertainty. 

The model uncertainties can be broadly grouped Into two categories: 

(1) uncertain external inputs 

- reference commands 

- environmental disturbances 

- biases or drifts in a failed sensor 

(2) uncertain internal dynamics 

- plant model errors 

- sensor failures 

- controller or estimator reconfigurations from DIA 

The representations of the model uncertainty vary according to how well 
its structure is known. For a highly structured representation, e.g., 
aerodynamic coefficients in flight control, engine model parameters in engine 
control, the uncertainty can be represented by defining a range of variation 
in the parameter space. For less structured uncertainties, bounds on errors 
can be defined. The model uncertainties can in general be classified as 
additive or multiplicative. The additive type of model uncertainty Is defined 
as follows 
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with 


P(s) = P Q (s, c *) + a(s) 


o[A(s)] < 6 (w) -V" w > 0 ( A. 1 ) 

a 

where P(s) is the transfer function of the plant, P (s, c*) Is the 

o 

nominal parameterized model of the plant, with structured uncertainty c*. In 

other words, P q (s, c *) 1s a known function of c*. but the values of c* 

are uncertain. The function a(s) is the variation in the nominal model, 

P , and is an unstructured uncertainty. It is unknown but limited to be 
o 

some function bounded by & (<j), where 6 (w) is a positive 

a a 

scalar function which defines a bound on P in the neighborhood of P. It 

o 

can be viewed as a frequency dependent "radius" of uncertainty of the true 
plant, P(s), about some model P q (s, c *) for c *- Figure A. la 
illustrates the additive type of uncertainty. In general a good model will be 
well known at low frequencies where 6 (w) will be small, and less well 

a 

known at high frequencies where 6 («) will be large. This type of 

a 

curve is characteri Stic of unmodeled, uncertain, high-frequency phenomena. 

Note that in equation (A.l), the structure of a(s) is not defined, and may 
be caused by a variety of mechanisms (for example, parameter changes, 
unmodeled dynamics, etc.). The two types of multiplicative uncertainties are 
the input multiplicative and the output multiplicative and are given by the 
following equations. 

Input multiplicative: 

P(s) = P q (s, c*) [I + A(s)J 

with 


<j[ A( s ) ] < <S (w) -V-u > 0 (A. 2) 

m 

Output multiplicative: 

P(s) = [I + Ms)] P q (s, c*) 

wi th 
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-v-w > 0 


(A. 3) 


a[a(s) ] < ^ m (<*>) 

Figure A. 1b and A.lc shows the representations of the multiplicative type of 
uncertainties . 

For the multivariable jet engine example [12], an output multiplicative 
bound was determined. This bound takes into account both the structured and 
unstructured uncertainties. This will represent a bound on such model 
uncertainties as, unmodeled dynamics, parameter variation in system matrices 
(A, B, C, and 0), unmodeled nonlinearities, reduced-order model ing, and 
linearization. 

A . 1 . 1 Why Determine Bounds on Model Errors? 

Under the NASA contracts NAS3-22481 and NAS3-23282, the feasibility of 
the OIA concept for application to the multi variable jet engine was 
demonstrated. However, several problem areas were identified which are stated 
below: 

(1) steady-state and dynamic mismatch of the simplified nonlinear 
models of the engine; 

(2) steady-state estimator errors, with no sensor failures; 

(3) instabilities when accommodating failures; 

(4) accommodation inaccuracies; and 

(5) missed detections and false alarms. 

These problems arise from system uncertainties. 

The fundamental control design problem is to establish the relationship 
between performance, robustness, and system uncertainty. The first two are 
conflicting requirements and need trade-off or design compromises to meet the 
system requi rements . A unified method of approach is to: 
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(1) establish quantitative relationships between performance, 
robustness, and system uncertainty; 

(2) recognize that the plant model error, sensor failures, and DIA 
reconfigurations all belong to the same class of system 
uncertainty. 

Figure A. 2 illustrates how the system uncertainties can be isolated. The 
dynamic uncertainties such as model errors, sensor failures, and DIA 
reconfigurations, are "inside" the system and function as a feedback loop 
around the "interconnection" system. This system maps the external and 
internal uncertainties into the outputs, i.e., tracking error and filter 
residuals. The dynamic uncertainties propagate in a specific way so as to 
cause a quantifiable uncertainty about the map from the input into the 
outputs, provided bounds can be found for the dynamic uncertainties. These 
bounds are obtainable from simple input/output system tests, and are to be 
used in robustifying the filter/s and determining the thresholds. 


A. 1.2 Model Error Bound for Output Multiplic ative Error 

A bound for model uncertainties can be obtained experimentally as shown 
in Figure A. 2. In Figure A. 2a, P represents a high-fidelity simulation of 
the plant whereas P Q is a simplified low-order linearized model of the 
plant. A sinusoidal°test input is applied to the plant and the model. The 

error is defined as 

e = P(s) Su - P (s, c*) Su 

= (I + A) P (s, c*) Su - P (s, c*) Su 
o o 

= AP (S, C*) SU (A - 4) 

0 

The normalized error, S(o>), provides a worst-case frequency-dependent 
bound, and is given by 

II e ii„ 

<$(o) = max 2 (A. 5) 

u II P Q (s, c*) Su n 2 
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Figure A. 2a Experiment to Determine Model Error Uncertainty 


*(") = IMI 'W y ml l wich u ■ V a si " 



1 10 100 
« FREQUENCY 


Figure A. 2b Frequency-Dependent Uncertainty 
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where IM^ is the norm, and 6(w) is the bound on model 
uncertainty. Figure A. 2b shows a typical bound on the model error 
uncertainties, which includes structured and unstructured uncertainties. 

A. 2 QUANTIFICATION OF MODEL ERROR UNCERTAINTIES 

This section describes a method to quantify the model error uncertainties 
using the experimental procedure given In Section A. 1.2. The error bound was 
determined for a jet engine model [12]. A description of this model Is also 
included in this section. 

A. 2.1 Computation of Model Error Bound using the Turbofan Model 

An experimental procedure to determine a bound on model error was 
discussed in Section A. 1.2. A nonlinear fourth-order engine model (HYTESS - A 
Hypothetical Turbofan Engine Simplified Simulation [12]) was used to represent 
the physical plant, l.e., the turbofan engine (see Figure A. 2a). The model 
P was represented by a linear model of the engine at sea level static 
flight condition and PLA at 36° (0/0/36). By applying the same sinusoidal 
input to P and P , an error signal between the outputs is computed using 
Eq. (A. 4). A frequency-dependent error bound is determined using Eq. (A. 5). 

The procedure for determining the error bounds Is shown in Figure A. 2. A 
step-by-step discussion follows. 

The states of the system are: 

xl = Fan Speed, SNFAN (Nl) - RPM 
x2 = Compressor Speed, SNCOM (N2) - RPM 
x3 = Burner Exit Slow Response Temperature, Tt41o - °R 
x4 = Fan Turbine Inlet Slow Response Temperature, Tt4.51o - °R 
The engine inputs are: 

ul = Main Burner Fuel Flow, WFMB - lb/hr 
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u2 = Nozzle Jet Area, Aj - ft? 

u3 = Fan Guide Vane Angle, FGV - deg 

u4 = Compressor Stator Vane Angle, SVA - deg 

u5 = Compressor Bleed Flow, BLC - % 

The engine outputs are: 

yl = Fan Speed, SNFAN (Nl) - RPM 

y2 = Compressor Speed, SNCOM (N2) - RPM 

y3 = Burner Pressure, PT4 - psia 

y4 = Augmentor Pressure, PT 6 - psia 

y5 = Fan Turbine Inlet Temperature, FTIT - °R 

The system matrices for the linear model and HYTESS are shown in Tables A.l 
and A. 2, respectively. 

A. 2. 2 Description of the Engine Models 

Two models of a jet engine are discussed in Section A. 2.1, namely, HYTESS 
and a linear model at 0/0/36. Both models are fourth-order state space 
models. HYTESS is a nonlinear model generated by scheduling linear model 
matrices (A, A 1 B, C and D) over the flight envelope using normalized 
variables (such as 6 , e, P/ 8 , T/e, N//e) as scheduling parameters. 

The linear model was extracted from HYTESS at the flight condition of altitude 
= 0 feet, Mach No. = 0, and PLA = 36°. 

The linear fourth-order state space models are of the form 
x = Ax -i- Bu 
y = Cx + Du 

where 

x = [ x ^ x 2 *2 * 4 ] 
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Table A.l 


Jet Engine System Matrices 


A 


-3. 9180 
-0. 2360 
-0. 0033 
-0. 0073 

2. 9550 
-2. 1480 
-0. 0043 
-0. 0146 

-1. 5450 
8. 4370 
-0. 6663 
-0. 0499 

6. 2260 
0. 0482 
0. 0003 
-2. 0000 


1. 0D+04 

* 




0. 0001 
0. 0001 
0. 0000 
0. 0000 

0 . 3071 
0 . 0395 
0 . 0002 
0 . 0004 

- 0 . 0068 
-0. 0003 
0. 0000 
0. 0000 

0. 0021 
0. 0015 
0. 0000 
0. 0000 

-1. 3290 
-0. 7886 
0. 0034 
0. 1618 

m 

1. 0000 
0. 0000 
0. 0147 
0. 0031 
-0. 0377 

0. 0000 
1. 0000 
o: 0284 
0. 0001 
-0. 0731 

0. 0000 
0. 0000 
-0. 0011 
-0. 0010 
-0. 2495 

0. 0000 
0. 0000 
-0. 0009 
-0. 0007 
0. 0018 


m 

0. 0000 
0. 0000 
0. 0087 
0. 0008 
0. 2050 

0. 0000 
0. 0000 
-8. 2350 
-5. 7450 
18. 5200 

0. 0000 
0. 0000 
0. 1824 
0. 0462 
-0. 5299 

0. 0000 
0. 0000 
0. 1855 
-0. 0032 
-1. 7080 

0. 0000 
0. 0000 
-175. 8000 
-6. 8890 
482. 9000 
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Table A. 2 


FJ. 1.-1) — *- 
F 12. 1 ) - 
F<3. 1 ) - 
F ( 4 , 1 ) - 
F < 1, 2>:..» 
F(2L21_- 
F(3.'2r*“ 
F(4,2> * 
F< 1. 3) - 
F(2, 3) « 
F(3. 3) - 
F < 4. 3) - 
F< 1,4) - 
F ( 2. 4) » 

F < 3. 4 ) “ 

F ( 4. 4 ) * ■ 


FICU, 1 ) 
FIC ( 2. 1 ) 
FIG<3. 1 ) 
FIG < 4. 1 ) 
FI G < 1. 2) 
FIG < 2» 2 ) 
FIC ( 3. 2) 
FIC < 4. 2) 
FIC< 1. 3) 
FIG ( 2. 3) 
FIC ( 3. 3) 
FIC ( 4. 3) 
FIC( I, 4) 
FIG < 2/ 4) 
FIC < 3» 4) 
FIC i 4. 4) 
FIC ( 1. 5) 
FIG < 2i 5) 
FIC<3. 5) 
FIC< 4. 3) 


Simplified Nonlinear Model Matrices 


FULL ENVELOPE MODEL OF A MATRIX 


^0. .9Ae3E^l^DL1^0.-1886E-2*XPT6**2-2. 463 
-0. 1 820E— 3*XN l +0: 27 68E-9*XPT 6* XN 1 **2+0. 7721 
-0. 2201 E— 2*XPT6+0. 8295E-4*XTT45-0. 8207E-1 
0. 6353E-2/DL1+0. 1O38E-0*XPT4**3-O. 2724E-1 
-O. 2596E+l*THl+0. 81 59E-1/DL1+5. 70 
— O.. 2381E+l/THlr:0_ 2348Et— 7*XN1**2+1.. 549 
-O. 4095 E— 2*DLrT-0^ 0 IT 4E-2/DL r+0. 121 2E- 1 
-0. 1010E-l*DLl-0. 2312E-1/DL1+0. 2996E-1 
-1 597 

-0. 1569/DL1+0. 23S2E— 1 *XPT4+4. 2890 
-0. 6377/DL1+0. 5218/TH1-0. 5229 
-0. 1156 

-0. 1 143/DL1+0. 8349E-7*XN1**2+1. 074 
0. 6013E-1 
0. 2632E-2 

-0. 1914E+1/DL1+0. 1564E+1/TH1-1. 568 


FULL ENVELOPE MODEL OF A -1 B MATRIX 


- 0. 9327E-3* XN 1 -0. 31 45E-1 1 *XN1 **3-6. 56 

- -0. 3616*TH1 -0. 3S50E+B/XN1 **2+0. 4541 

0. 1 1 83E-3*XPT4-0. 363E-9*XPT4**3-0 372E-1 
0. 9333E-3*XPT6-0. 9744E-9*XNl*XPT6**2-0. 306E-1 
-0. 2971E+2*XPT6+0. 21B2E+2/ i TH1*DL1 >-260. 4 
0. 1121 *XNl-0. 1 545E-2*XPT4«XPT6**2-721 . 8 
0. 584 5E— 1+XPT4— 0. 7198E-4*XPT4*XPT6**2-1 568 
0. 3761 *XPT6— O. 5675E— 4*XPT4*XPT6**2— 1'. 418 
0. 1627*XPT4— O. 4 133E— 4*XPT4*XPT6**2— 2. 104 
0. 275*XPT6-6. 304 
0. 1306E— 1*XPT6— 0. 2379 
0. 1068E-1 *XPT6-0. 1881 
0. S670E+1/TH1+0. 5090E-l*XPT4-26. 19 
0. 1 16SE+3*THl+0. 1 092E+3/TH1 — 220 7 
0. 2010*TH1-0. l262E-2*XPT4+0. 2459 
0. 9386E- 1 *TH 1 — O. 1 1 57E-2*XP74+C. 2951 
0. 5357E+4/ ( DL1 *RTH1 ) -75. 87 
0. 3754E+4/ < DL1 *RTH1 > + 157. 9 
-O. 9392E+2*<RTH1/DL1 )— O. 7399E+1 *» THl /EL 1 **2 ) +8 
-0. 8524E+2*(RTH1/DL1 ) -0. 3902E+1 *< TH1 /DL 1 **2 1 +9. 
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Table A. 2 (Continued) 


ORIGINAL PASS u> 
0Z POOR QUALITY 


FULL ENVELOPE MODEL OF C MATRIX 


HI I. I I - «. 

h< 2 . i » - o. 

HO. || - -O. 281 46*1 /KPT4-0. 23786-3*lPT4.«»T4*0. *3416-1 

Ml 4* 1) - -O. 3DI8/XPT4-0. 19306-1 0*IH1 **2*0. **346-2 

MIX II - 0. 87146-3**PT«*0 I704£*2/IPT4-O. 3047 

MIX 1 1 - O. 91B4/TH1-0. 349|C*4/rPT4**2*0. 1318 

Ml 7. II - O. 31l9C-2/*PT4-0. 29796-t /lPT***2-». *26*6-4 

HI 4*3) - O. 

HO. 21 - i. 

HO* 2) 0..291 3€r l/TH! *0. 44446-3* *PT*-0. 12336-1 

11(1.31*-'^. 90006- 7* *M 1-0. 4X446-1 /IPT4-0. 1 1 326-2 

HO* 21 - -O. l3276-3**TT43-0. 19076*4/1M1*C. *014 

HI 4. 21 - O. 32876-1 

Ml 7. 21 - O. 2«99€-3/0Cl*0. 4102/K*T4~2~0. 3*236—4 
HI I • 3 1 • 0. 

HO. 3) - 0. 

Hi 3* 3 1 - -0. 73706-7* **T4* *3-0. 48476-3 

HI 4. 31 - -0. 32046-4 *XN2*0. 23996-2 

HO* 31 - 0. 31116-1*1X1*0. 34846 -4*UT43-0. 3*12 

MIX 3) - 0. 14436-i/OCi-O. I2716-3*IHI*0. 3134 

Ml 7. 31 - -*0. 13226-3 

MU. 4) - 0. 

MO. 41 - 0. 

MIX 41 - -0. 1 9806-3* JCPT4*0. 31736-2 
MI4. 4) - -0. 43006-3* *HT4-0. 310 26-4 
MIX 4) - 0. 84846-2 

HI 4. 41 - 0. 11406* I /DC 1-0. 13006- 3* CM 1*0. 31*3 

HIT. 41 - -0. 13246*3 


FULL ENVELOPE MODEL OF D MATRIX 


Dll. 1 1 - O. 

0 ( 2 . 1 1 - o. _ 

OIXW - - 0 . 41336-14*KW1**3*0. l983£-J*Cpr4**2*<J *9726-2 
014. II - -O. 94426-7*CNi-0. 22076*3/7*1 « *2*0. 20^46-2 

01X41 - -O. 334 16-1*1X1*0. 3133€*3/X*T4*0. 49*46-1 

014* II - -0. 48 306*4 /IN 1 *~0. 1 32 46 *-?/KK2* * 2*0. 3303 
017. II - -O. 83296-2/7PT4*0. 10336-4 
011.21 - 0. 


01X2) - 0. 

01X21 - -0. 7192£-3*IPT4**3*1. 233 

014.21 - -0. 4047*CPT4*2. 900 

00*21 - 0. 2438€-2*IPTA**3-7. U9 

014*21 - -0. 27046*4/1X1*0. 83346*7/1X1*131. 4 

017.21 - -0. 39426-4* CPT4**3*0. 20296-3 

011.31 - 0. 

01X31 - 0. 

01X31 - 0. 3334C-2*KPT4-0. 1470€-3*f^T4»lPT4**2-0. 4 401 

014. 31 - 0. 881 36-3* 7PT4-0. 28346-8*CFT4«O-0. 44386-1 

01X31 - 0. It l0C-4*t^T4**3-0. 17«9«-3*INI**PT4*0. 3499 

014*31 - 0. 13476*1 *KPT4- 1 2. 47 

017.31 - 0.20496-3 

011.41 - 0. 

OCX 41 - 0. 

01X41 - 0. 7838*1X1-0.4444 

014.41 - -0.34306-2 

01X41 - -0. 21476*t*TXt*0. 271 AC-7*fMl**2-t. 34 
01X41 - 1.073 

017.41 - 0. 33276-2* TX 1-0. 84346-2 

DU. 31 - 0. 

OCX 3) - 0. 

01X31 - -0. 31 43€*3/0C 1-0. 1 4836*1 74*423. 9 

OCX 31 - -0. |444C*2*<*THt/DLl 1-0. 93486-1* ***4 *2 4. 09 

01X31 - 0. 11036*4/011-431. 2 

01X3) • -O. 74826 *4 /0L1* 1434. 0 

017.31 * 0. 9834/00*0. 33806*1 
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U T = [u 1 u 2 u 3 U 4 ] 

y T = [y-, y 2 v 3 v 4 ] 

A. 3 DETERMINATION OF MODEL ERROR BOUND 

The procedure for computation of the error bound on modeling error is 
shown in Figure A. 3, and a step by step discussion follows. 

Step 1 : The Input to the system is a sinusoidal signal superimposed on a 

constant signal. This input Is of a specific frequency and of high amplitude, 
i.e., the amplitude of the sinusoid should be carefully selected so that the 
region of operation of the simplified model remains linear. The amplitudes of 
the sinusoid for each input were determined from Ref. 24. These input 
amplitudes give the best overall match between the linear and the nonlinear 
models of the engine. These amplitudes are as follows: 


u(l) 

= WFMBH 

3 

% 

u(2) 

= AJ 

3 

% 

u( 3) 

= CIVV 

5 

degrees 

U(4) 

= RCVV 

1 

degree 

u(5) 

= BLC 

0 

.2 % 


The inputs were excited one at a time according to the procedure. The outputs 

of P and P are sinusoids superimposed over a transient (system 
0 

transient) , as shown in Figure A. 4. Once the transient dies out, the outputs 
are sinusoids of constant amplitude. The magnitude of the outputs are 
determined by computing the amplitude of the sine waves. The error, e, is 
determined by subtracting the output of the plant, y, from the output of the 
model, ym, and then computing the magnitude. 

Step 2 : The bound, 4(o>), on the model error is given by 
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Figure A. 3 Procedure for Computing Error Bound, 6(«) 
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M«) = 


(A. 6) 


4 n 4 i5 

4 51 4 55 

Each input when excited, produces errors in all the five outputs, which 
determines one column of the 4 matrix. For example, if u^ is 
modulated with a sinusoid, it produces 

4-j-! - ny-j - ynyi / Iiynyi, « 21 = Iiy 2 - ynyi / Iiym 2 ii, ... 

4 = lly 5 " yf V 7 ,ly V- 

This gives the first column of the 6 matrix. Similarly, other four inputs 

are excited one at a time to determine the other four columns of the 4 
matrix. To determine the bound at a given frequency, 25 elements of the 4 
matrix have to be computed. 

Step 3 : The simplified model has a bandwidth of about 5 hz. A range of 

frequencies from 0 hz (OC) to 9 hz was chosen for computation of the 4 
matrix. The discrete frequencies at which 4 was computed are 0.1, 0.3, 
0.5, 0.7, 0.8, 1., 2., 3., 4., 5., 7., 9. hz. 

Step 4: At each frequency, the maximum singular value of the 4(w) 

matrix is computed. This 1$ denoted by o[4(<*>)], and is the worst case 

bound at that frequency. Figure A. 4 illustrates a plot of a against 
frequency. 

A. 4 DISCUSSION 

The bound on the model errors determined in this procedure has some 
limitations. The plant Is represented by a simplified nonlinear simulation 
and the model is linearized at one flight point. This limits the validity of 
the model error bound to only one flight condition. This also does not 
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Include the errors which arise from the unmodeled nonlinearities and the 
linearization of P(s) to produce simplified model P Q ( s ) . For a more 
general error bound, the plant should be represented by the cycle deck, and 
the bound 5 should be determined at all the representative flight points 
in the flight envelope. However, cost is an important consideration in such 
an exercise. 

There are marked differences between the theoretical curves of Figure A. 2 
and Figure A. 5. The error bound Is expected to grow at high frequencies 
(Figure A. 2b), but as seen in Figure A. 4, it is constant. This Is attributed 
to the fact that the plant and the model are of the same order and have 
approximately the same A, B, C, and D matrices at 0/0/36. The difference in 
the system matrices is that the model is represented by the constant A, B, C, 
and D matrices (Table A.l), whereas the plant P Q (s, c *) represented by 
the A, B, C, and 0 matrices which are polynomial functions. The error bound 
in this case is a constant. This can be shown using equation A. 3 as follows, 

P(s) = [I + A( s) ] P (s, c*) (A. 7) 

o 

where P(s) and P (s, c*) are of the same order. In reality, a(s) 

0 

would be a complicated transfer function of order different from the order of 
P(s) or P o (s, c*). In this case it is simply a constant and therefore the 

error bound a(u>) is a constant in Figure A. 5 as shown by the following 
equation, 

nP(s) - P (s, C*)ll 

lim rjrfrnj = ha(s)ii - constant (A. 8) 

s -* ® o' 5 ' 

It is expected that if high-frequency dynamics is added to the plant P(s), the 

high frequency response of a in Figure A. 5 will follow the pattern of 
Figure A. 2b. 

The difference between the two models P(s) and P Q (s, c *)» causes the 
peak at low frequency in Figure A. 5. This difference gives rise to different 
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cf[ 6(a)) 



3 . 4 . 5 . 6 . 7 . 8 . 9 

f (hz) 


Figure A. 5 Maximum Singular Value Plotted Against Frequency 
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dynamic responses. It Is expected that if the two models are exactly the same 

(no modeling errors), the plot of a should be a constant, l.e., o(w) ■ a 
(constant) = 0. 

A. 5 SUMMARY 

A procedure to compute the bound on the modeling errors is developed and 
demonstrated on an engine example. The bound is limited to only one flight 
condition. For a more general bound, a comprehensive experiment will have to 
be run where the plant is represented by the cycle deck and the experiment 
conducted at a number of operating points on the flight envelope, with 
particular emphasis on the choice of the amplitude of the input signals. 
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